Brain Age Prediction code
Dataset
- Training HC (N=541)
- KU HC (N=209)
- CNUH patient (N=233)
- TRS : 46
- FE-SSD : 72
- Input data : the 77 gray matter FreeSurfer features.
- ICV : 1 variable
- Subcortical volume : 8 variables
- Mean cortical thickness : 34 variables
- Mean cortical surface area : 34 variables
Total 77 variables
- ICV : 1 variable
Strategy of brain age prediction
Sample exclusion criteria
- Exclude samples below 18 years and above 75 years old.
- Exclude samples included in scanning sites with less than 20 control samples.
Data partitioning
- Train HC Control : Healthy controls from CNUH+Beijing.
- Test HC Control : Healthy control from KU.
- Test Patient : Patient from CNUH.
Model development
- Perform Combat to remove site effect.
- After removing site effect, build brain age prediction model based on the 77 gray matter FreeSurfer features using Ridge, SVR, RVR, 10-fold Nested cross-validation (10 time).
Test prediction model in Test control data, CNUH patient using variance explained on cross-validation (VEcv), Mean absoulte error (MAE), Peason correlation coefficient (\(r\)).
Calculate Brain-PAD (predicted Age - chronological Age) and corrected Brain-PAD (Goodness-of-fit test)
Data Process
library(openxlsx)
library(dplyr)##
## 다음의 패키지를 부착합니다: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionwhole_dat <- read.xlsx("2023-1114-BrainAge-sMRI-Final.xlsx",sheet=1,startRow = 1)Sample exclusion with Age & Site
# Age (filtering subject with age <18 or >75)
whole_dat$Age <- as.numeric(whole_dat$Age)
whole_dat$Age <- floor(whole_dat$Age)
sum(whole_dat$Age < 18 | whole_dat$Age > 75)## [1] 93table(whole_dat$Site[which(whole_dat$Age < 18 | whole_dat$Age > 75)])##
## IPCAS SWU
## 92 1whole_dat <- whole_dat[-which(whole_dat$Age < 18 | whole_dat$Age > 75),] #1075
# Site
table(whole_dat$Site)##
## BNU CNUH HNU IACAS IPCAS JHNU KU SWU XHCUMS
## 56 523 3 2 100 30 214 70 23whole_dat <- whole_dat[-which(whole_dat$Site=="HNU" | whole_dat$Site=="IACAS"),]
table(whole_dat$Site_2)##
## BEI CNUH KU
## 279 523 214table(whole_dat$Status)##
## Control Patient
## 783 233# Education
whole_dat$Edu[whole_dat$Edu==1 & !is.na(whole_dat$Edu)] <- 7.5
whole_dat$Edu[whole_dat$Edu==2 & !is.na(whole_dat$Edu)] <- 14
whole_dat$Edu[whole_dat$Edu==3 & !is.na(whole_dat$Edu)] <- 18
whole_dat$Edu <- as.numeric(whole_dat$Edu)Reduce and seperation Reduce of dataset(KU reslut reference)
table(whole_dat$Status)##
## Control Patient
## 783 233table(whole_dat$Site_2)##
## BEI CNUH KU
## 279 523 214'%!in%' <- Negate('%in%')
whole_dat <- whole_dat %>% filter(whole_dat$Freesurfer_Number %!in% c("H148","H252","J373"))
whole_dat <- whole_dat %>% filter(whole_dat$Freesurfer_Number %!in% c("K33","K49","K65","K75","K216"))
HC_data <- whole_dat[whole_dat$Status=="Control",]
Patient_data <- setdiff(whole_dat,HC_data)
HC_data <- HC_data[,c(-2)]
Patient_data <- Patient_data[,c(-2)]Combat
library(sva)
library(factoextra)
whole_dat$Site <- ifelse(whole_dat$Site=="CNUH","JBUH",whole_dat$Site)
whole_dat$Site <- ifelse(whole_dat$Site=="KU","KUAH",whole_dat$Site)
batch <- as.numeric(as.factor(whole_dat$Site))
mod <- model.matrix(~1+Age+Sex+Status,data=whole_dat)
Combat_dat <- ComBat(dat=t(whole_dat[,9:85]), batch=batch, mod=mod)
Combat_dat <- t(Combat_dat)
Combat_dat <- data.frame(SubjID=whole_dat$Freesurfer_Number,Age=whole_dat$Age,
Sex=whole_dat$Sex,Site=whole_dat$Site,Combat_dat)
Train_HC_data <- Combat_dat[grep("H|J",whole_dat$Freesurfer_Number),]
Test_HC_data <- Combat_dat[grep("K",whole_dat$Freesurfer_Number),]
Test_Patient_dat <- Combat_dat[grep("C|T",whole_dat$Freesurfer_Number),]
Patient_clinical <- read.xlsx("Patient_clinical_292_0630.xlsx",startRow=1)
Patient_clinical <- Patient_clinical[match(Test_Patient_dat$SubjID,Patient_clinical$freesurfer_ID),]
ID_FESSD <- Patient_clinical$freesurfer_ID[!is.na(Patient_clinical$`FE-SSD`)]
ID_schizo <- Patient_clinical$freesurfer_ID[Patient_clinical$Diagnosis=="schizophrenia"]
Test_TRS_dat <- Test_Patient_dat[grep("T",Test_Patient_dat$SubjID),]
Test_FESSD_dat <- Test_Patient_dat[Test_Patient_dat$SubjID %in% ID_FESSD,]
Test_schizophrenia_dat <- Test_Patient_dat[Test_Patient_dat$SubjID %in% ID_schizo,]
par(mfrow=c(3,2))
plot(density(Train_HC_data$Age,from=15,to=80),main="Age of Train control data")
plot(density(Test_HC_data$Age,from=15,to=80),main="Age of Test control data")
plot(density(Test_Patient_dat$Age,from=15,to=80),main="Age of Test patient data")
plot(density(Test_schizophrenia_dat$Age,from=15,to=80),main="Age of Schizophrenia data")
plot(density(Test_FESSD_dat$Age,from=15,to=80),main="Age of Test FE-SSD data")
plot(density(Test_TRS_dat$Age,from=15,to=80),main="Age of Test TRS data")
PC <- prcomp(scale(whole_dat[,9:85]))
A <- fviz_pca_ind(PC, label="none", habillage=whole_dat$Site,
addEllipses=TRUE, ellipse.level=0.95, palette = "Dark2",title="Before adjustment") +
labs(x="Dimension1 (29.0%)", y="Dimension2 (20.7%)")
PC <- prcomp(scale(Combat_dat[,5:81]))
B <- fviz_pca_ind(PC, label="none", habillage=Combat_dat$Site,
addEllipses=TRUE, ellipse.level=0.95, palette = "Dark2",title="After adjustment") +
labs(x="Dimension1 (30.0%)", y="Dimension2 (20.8%)")
library(MASS)
f <- as.formula(paste("Site~",paste(colnames(whole_dat)[9:85],collapse ="+"),sep=""))
lda <- lda(f, data = whole_dat)
lda.value <- predict(lda)
lda.data <- data.frame(Site = whole_dat$Site, lda = lda.value$x)
C <- ggplot(lda.data) + geom_point(aes(lda.LD1, lda.LD2, colour = Site), size = 2.5)
lda <- lda(f, data = Combat_dat)
lda.value <- predict(lda)
lda.data <- data.frame(Site = Combat_dat$Site, lda = lda.value$x)
D <- ggplot(lda.data) + geom_point(aes(lda.LD1, lda.LD2, colour = Site), size = 2.5)
A;B;C;D
library(ggpubr)
png("s1.png", width = 1200, height = 800)
ggarrange(A,B,C,D,
labels = c("A","B","C","D"),
ncol = 2, nrow = 2,
font.label = list(size=20))
dev.off()## png
## 2# PCA outlier dection in Train HC data
PC <- prcomp(scale(Train_HC_data[,5:81]))
fviz_pca_ind(PC,col.ind = "#00AFBB",repel=TRUE)## Warning: ggrepel: 465 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
fviz_pca_ind(PC, addEllipses=TRUE, ellipse.level=0.95, palette = "Dark2", title="") + labs(x="Dimension1 (30.3%)", y="Dimension2 (18.1%)")
out <- c(631,95,171,75,233,646,159,766,223,247,267,5,154,186,192,268,36,238,256,3,249,41,211,197,657) #1067
Train_HC_data$SubjID[rownames(Train_HC_data) %in% out]## [1] "H03" "H05" "H36" "H41" "H75" "H95" "H155" "H160" "H172" "H187"
## [11] "H193" "H198" "H212" "H224" "H234" "H239" "H248" "H268" "H250" "H260"
## [21] "H263" "J161" "J248" "J259" "J375"length(out) # 25## [1] 25sum(!rownames(Train_HC_data) %in% out) # 566 --> 541## [1] 541Train_HC_data <- Train_HC_data[which(!rownames(Train_HC_data) %in% out),]Demographics and Clinical characteristics
#Train_HC vs Test_HC vs Schizophrenia
summary_dat_3 <- data.frame(rbind(Train_HC_data,Test_HC_data,Test_schizophrenia_dat))
Edu_dat <- whole_dat[,c(1,6)]
colnames(Edu_dat)[1] <- c("SubjID")
summary_dat_3 <- inner_join(summary_dat_3,Edu_dat,by="SubjID")
summary_dat_3$set <- rep(c(1:3),c(nrow(Train_HC_data),nrow(Test_HC_data),nrow(Test_schizophrenia_dat)))
summary_dat_3$set <- as.factor(summary_dat_3$set)
three.fun_var <- function(var,group,DATA) {
a11 <- paste("n=",sum(!is.na(DATA[var][DATA[group]==1])),sep="")
a21 <- paste(round(mean(DATA[var][DATA[group]==1],na.rm = T),2),
round(sd(DATA[var][DATA[group]==1],na.rm = T),2),sep = " ± ")
a_1 <- paste(a21,"(",a11,")",sep="")
a12 <- paste("n=",sum(!is.na(DATA[var][DATA[group]==2])),sep="")
a22 <- paste(round(mean(DATA[var][DATA[group]==2],na.rm = T),2),
round(sd(DATA[var][DATA[group]==2],na.rm = T),2),sep = " ± ")
a_2 <- paste(a22,"(",a12,")",sep="")
a13 <- paste("n=",sum(!is.na(DATA[var][DATA[group]==3])),sep="")
a23 <- paste(round(mean(DATA[var][DATA[group]==3],na.rm = T),2),
round(sd(DATA[var][DATA[group]==3],na.rm = T),2),sep = " ± ")
a_3 <- paste(a23,"(",a13,")",sep="")
f <- as.formula(paste(var,"~",group))
Aov <- aov(lm(f,data=DATA))
b <- round(summary(Aov)[[1]][["Pr(>F)"]][1],4)
y <- DATA[var][,1]
group <-DATA[group][,1]
posthoc <- pairwise.t.test(y, group, "bonferroni")
posthoc <- data.frame(posthoc[3])
c <- round(posthoc[1,1],4)
d <- round(posthoc[2,1],4)
e <- round(posthoc[2,2],4)
return(matrix(c(a_1,a_2,a_3,b,c,d,e),1,7))
}
TABLE_clinical_3 <- three.fun_var("Age","set",summary_dat_3)
TABLE_clinical_3 <- rbind(TABLE_clinical_3,three.fun_var("Edu","set",summary_dat_3))
t <- table(summary_dat_3$Sex,summary_dat_3$set)
Sex <- c(paste(t[,1],collapse = "/"),paste(t[,2],collapse = "/"),paste(t[,3],collapse = "/"),round(chisq.test(t)$p.value,4))
Sex <- c(Sex,NA,NA,NA)
TABLE_clinical_3 <- rbind(Sex,TABLE_clinical_3)
rownames(TABLE_clinical_3) <- c("Sex(M/F)","Age","Education, years")
colnames(TABLE_clinical_3) <- c("Health control (CNUH, Beijing)","Health control (KU)","Schizophrenia (CNUH)","p-value","Group 1 vs Group 2", "Group 1 vs Group 3","Group 2 vs Group 3")
TABLE_clinical_3## Health control (CNUH, Beijing) Health control (KU)
## Sex(M/F) "281/260" "77/132"
## Age "31.9 ± 13.7(n=541)" "37.66 ± 14.66(n=209)"
## Education, years "13.94 ± 2.41(n=267)" "13.83 ± 2.35(n=209)"
## Schizophrenia (CNUH) p-value Group 1 vs Group 2
## Sex(M/F) "108/88" "2e-04" NA
## Age "39.34 ± 12.9(n=196)" "0" "0"
## Education, years "13.19 ± 3.17(n=195)" "0.0074" "1"
## Group 1 vs Group 3 Group 2 vs Group 3
## Sex(M/F) NA NA
## Age "0" "0.6544"
## Education, years "0.0084" "0.0473"#TestHC vs FE-SSD vs TRS
summary_dat_p3 <- data.frame(rbind(Test_HC_data,Test_FESSD_dat,Test_TRS_dat))
summary_dat_p3 <- inner_join(summary_dat_p3,Edu_dat,by="SubjID")
summary_dat_p3$set <- rep(c(1:3),c(nrow(Test_HC_data),nrow(Test_FESSD_dat),nrow(Test_TRS_dat)))
summary_dat_p3$set <- as.factor(summary_dat_p3$set)
TABLE_clinical_p3 <- three.fun_var("Age","set",summary_dat_p3)
TABLE_clinical_p3 <- rbind(TABLE_clinical_p3,three.fun_var("Edu","set",summary_dat_p3))
t <- table(summary_dat_p3$Sex,summary_dat_p3$set)
Sex <- c(paste(t[,1],collapse = "/"),paste(t[,2],collapse = "/"),paste(t[,3],collapse = "/"),round(chisq.test(t)$p.value,3))
Sex <- c(Sex,NA,NA,NA)
TABLE_clinical_p3 <- rbind(Sex,TABLE_clinical_p3)
rownames(TABLE_clinical_p3) <- c("Sex(M/F)","Age","Education, years")
colnames(TABLE_clinical_p3) <- c("Test HC","FE-SSD","TRS","p-value","Group 1 vs Group 2", "Group 1 vs Group 3","Group 2 vs Group 3")
TABLE_clinical_p3## Test HC FE-SSD
## Sex(M/F) "77/132" "28/44"
## Age "37.66 ± 14.66(n=209)" "29.61 ± 12.26(n=72)"
## Education, years "13.83 ± 2.35(n=209)" "13.26 ± 2.79(n=72)"
## TRS p-value Group 1 vs Group 2
## Sex(M/F) "30/16" "0.002" NA
## Age "42.61 ± 9.9(n=46)" "0" "1e-04"
## Education, years "13.73 ± 2.36(n=46)" "0.2405" "0.277"
## Group 1 vs Group 3 Group 2 vs Group 3
## Sex(M/F) NA NA
## Age "0.0776" "0"
## Education, years "1" "0.9508"Comparison of the 77 gray matter Freesurfer features
three.fun_var <- function(var,group,DATA) {
a11 <- paste("n=",sum(!is.na(DATA[var][DATA[group]==1])),sep="")
a21 <- paste(round(mean(DATA[var][DATA[group]==1],na.rm = T),2),
round(sd(DATA[var][DATA[group]==1],na.rm = T),2),sep = " ± ")
a_1 <- paste(a21,"(",a11,")",sep="")
a12 <- paste("n=",sum(!is.na(DATA[var][DATA[group]==2])),sep="")
a22 <- paste(round(mean(DATA[var][DATA[group]==2],na.rm = T),2),
round(sd(DATA[var][DATA[group]==2],na.rm = T),2),sep = " ± ")
a_2 <- paste(a22,"(",a12,")",sep="")
a13 <- paste("n=",sum(!is.na(DATA[var][DATA[group]==3])),sep="")
a23 <- paste(round(mean(DATA[var][DATA[group]==3],na.rm = T),2),
round(sd(DATA[var][DATA[group]==3],na.rm = T),2),sep = " ± ")
a_3 <- paste(a23,"(",a13,")",sep="")
f <- as.formula(paste(var,"~","set"))
Aov <- aov(lm(f,data=DATA))
b <- round(summary(Aov)[[1]][["Pr(>F)"]][1],4)
y <- DATA[var][,1]
group <-DATA[group][,1]
posthoc <- pairwise.t.test(y, group, "bonferroni")
posthoc <- data.frame(posthoc[3])
c <- round(posthoc[1,1],4)
d <- round(posthoc[2,1],4)
e <- round(posthoc[2,2],4)
f1 <- as.formula(paste(var,"~","Age + Sex + Edu + set"))
Aov1 <- aov(lm(f1,data=DATA))
f <- round(summary(Aov1)[[1]][["Pr(>F)"]][4],4)
return(matrix(c(a_1,a_2,a_3,b,c,d,e,f),1,8))
}
conti_var <- setdiff(colnames(summary_dat_3),c("SubjID","Age","Site","Sex","set","Edu"))
TABLE_sMRI_3group <- NULL
for (i in 1:length(conti_var)){
TABLE_sMRI_3group <- rbind(TABLE_sMRI_3group,three.fun_var(conti_var[i],"set",summary_dat_3))
}
rownames(TABLE_sMRI_3group) <- conti_var
colnames(TABLE_sMRI_3group) <- c("Health control (CNUH, Beijing)","Health control (KU)","Schizophrenia (CNUH)","p-value","Group 1 vs Group 2", "Group 1 vs Group 3","Group 2 vs Group 3","age.sex.edu.adj-p")
TABLE_sMRI_3group <- data.frame(TABLE_sMRI_3group)
TABLE_sMRI_3group$p.value <- as.numeric(paste(TABLE_sMRI_3group$p.value))
TABLE_sMRI_3group$bonf.adj.p.value <- p.adjust(TABLE_sMRI_3group$p.value,method="bonferroni")
TABLE_sMRI_3group$bonf.age.sex.edu.adj.p.value <- p.adjust(TABLE_sMRI_3group$age.sex.edu.adj.p,method="bonferroni")
TABLE_sMRI_3group## Health.control..CNUH..Beijing.
## Avg_bankssts_thickavg 2.62 ± 0.13(n=541)
## Avg_caudalanteriorcingulate_thickavg 2.71 ± 0.18(n=541)
## Avg_caudalmiddlefrontal_thickavg 2.66 ± 0.11(n=541)
## Avg_cuneus_thickavg 1.9 ± 0.11(n=541)
## Avg_entorhinal_thickavg 3.47 ± 0.25(n=541)
## Avg_fusiform_thickavg 2.8 ± 0.11(n=541)
## Avg_inferiorparietal_thickavg 2.57 ± 0.11(n=541)
## Avg_inferiortemporal_thickavg 2.87 ± 0.11(n=541)
## Avg_isthmuscingulate_thickavg 2.39 ± 0.16(n=541)
## Avg_lateraloccipital_thickavg 2.3 ± 0.1(n=541)
## Avg_lateralorbitofrontal_thickavg 2.74 ± 0.12(n=541)
## Avg_lingual_thickavg 2.01 ± 0.11(n=541)
## Avg_medialorbitofrontal_thickavg 2.58 ± 0.13(n=541)
## Avg_middletemporal_thickavg 2.97 ± 0.12(n=541)
## Avg_parahippocampal_thickavg 2.69 ± 0.21(n=541)
## Avg_paracentral_thickavg 2.51 ± 0.13(n=541)
## Avg_parsopercularis_thickavg 2.68 ± 0.13(n=541)
## Avg_parsorbitalis_thickavg 2.8 ± 0.15(n=541)
## Avg_parstriangularis_thickavg 2.58 ± 0.12(n=541)
## Avg_pericalcarine_thickavg 1.61 ± 0.11(n=541)
## Avg_postcentral_thickavg 2.12 ± 0.11(n=541)
## Avg_posteriorcingulate_thickavg 2.57 ± 0.12(n=541)
## Avg_precentral_thickavg 2.62 ± 0.12(n=541)
## Avg_precuneus_thickavg 2.45 ± 0.11(n=541)
## Avg_rostralanteriorcingulate_thickavg 2.95 ± 0.16(n=541)
## Avg_rostralmiddlefrontal_thickavg 2.5 ± 0.11(n=541)
## Avg_superiorfrontal_thickavg 2.88 ± 0.12(n=541)
## Avg_superiorparietal_thickavg 2.28 ± 0.1(n=541)
## Avg_superiortemporal_thickavg 2.88 ± 0.14(n=541)
## Avg_supramarginal_thickavg 2.62 ± 0.12(n=541)
## Avg_frontalpole_thickavg 2.97 ± 0.24(n=541)
## Avg_temporalpole_thickavg 3.75 ± 0.22(n=541)
## Avg_transversetemporal_thickavg 2.34 ± 0.17(n=541)
## Avg_insula_thickavg 3.04 ± 0.14(n=541)
## Avg_bankssts_surfavg 967.57 ± 121.79(n=541)
## Avg_caudalanteriorcingulate_surfavg 645.48 ± 108.22(n=541)
## Avg_caudalmiddlefrontal_surfavg 2165.97 ± 304.8(n=541)
## Avg_cuneus_surfavg 1625.94 ± 191.04(n=541)
## Avg_entorhinal_surfavg 419.46 ± 68.21(n=541)
## Avg_fusiform_surfavg 3129.33 ± 324.28(n=541)
## Avg_inferiorparietal_surfavg 5216.56 ± 635.7(n=541)
## Avg_inferiortemporal_surfavg 3494.81 ± 429.24(n=541)
## Avg_isthmuscingulate_surfavg 1008.34 ± 133.08(n=541)
## Avg_lateraloccipital_surfavg 5117.72 ± 579.53(n=541)
## Avg_lateralorbitofrontal_surfavg 2588.52 ± 250.91(n=541)
## Avg_lingual_surfavg 3231.03 ± 386.93(n=541)
## Avg_medialorbitofrontal_surfavg 1914.23 ± 188.57(n=541)
## Avg_middletemporal_surfavg 3476.88 ± 401.7(n=541)
## Avg_parahippocampal_surfavg 662.46 ± 69.9(n=541)
## Avg_paracentral_surfavg 1432.7 ± 147.16(n=541)
## Avg_parsopercularis_surfavg 1504.98 ± 196.49(n=541)
## Avg_parsorbitalis_surfavg 728.13 ± 80.8(n=541)
## Avg_parstriangularis_surfavg 1450.58 ± 191.35(n=541)
## Avg_pericalcarine_surfavg 1565.86 ± 244.81(n=541)
## Avg_postcentral_surfavg 4197.01 ± 426.36(n=541)
## Avg_posteriorcingulate_surfavg 1202.75 ± 150.23(n=541)
## Avg_precentral_surfavg 4980.56 ± 461.28(n=541)
## Avg_precuneus_surfavg 4073.36 ± 461.12(n=541)
## Avg_rostralanteriorcingulate_surfavg 707.86 ± 109.49(n=541)
## Avg_rostralmiddlefrontal_surfavg 5662.36 ± 719.79(n=541)
## Avg_superiorfrontal_surfavg 6853.75 ± 716.93(n=541)
## Avg_superiorparietal_surfavg 5494.47 ± 593.85(n=541)
## Avg_superiortemporal_surfavg 3877.44 ± 387.67(n=541)
## Avg_supramarginal_surfavg 3952.15 ± 525.93(n=541)
## Avg_frontalpole_surfavg 277.22 ± 26.03(n=541)
## Avg_temporalpole_surfavg 473.96 ± 49.44(n=541)
## Avg_transversetemporal_surfavg 412.78 ± 52.03(n=541)
## Avg_insula_surfavg 2458.61 ± 232.1(n=541)
## ICV 1510324.43 ± 157050.4(n=541)
## Avg_LatVent 7516.06 ± 4127.39(n=541)
## Avg_thal 8073.94 ± 842.81(n=541)
## Avg_caud 3554.28 ± 456.93(n=541)
## Avg_put 5148.74 ± 589.6(n=541)
## Avg_pal 2020.69 ± 230.34(n=541)
## Avg_hippo 4268.04 ± 370.26(n=541)
## Avg_amyg 1758.63 ± 188.75(n=541)
## Avg_accumb 484.49 ± 81.47(n=541)
## Health.control..KU.
## Avg_bankssts_thickavg 2.6 ± 0.13(n=209)
## Avg_caudalanteriorcingulate_thickavg 2.68 ± 0.18(n=209)
## Avg_caudalmiddlefrontal_thickavg 2.64 ± 0.13(n=209)
## Avg_cuneus_thickavg 1.88 ± 0.11(n=209)
## Avg_entorhinal_thickavg 3.46 ± 0.26(n=209)
## Avg_fusiform_thickavg 2.78 ± 0.13(n=209)
## Avg_inferiorparietal_thickavg 2.55 ± 0.11(n=209)
## Avg_inferiortemporal_thickavg 2.85 ± 0.12(n=209)
## Avg_isthmuscingulate_thickavg 2.36 ± 0.17(n=209)
## Avg_lateraloccipital_thickavg 2.29 ± 0.11(n=209)
## Avg_lateralorbitofrontal_thickavg 2.71 ± 0.15(n=209)
## Avg_lingual_thickavg 1.99 ± 0.11(n=209)
## Avg_medialorbitofrontal_thickavg 2.55 ± 0.15(n=209)
## Avg_middletemporal_thickavg 2.95 ± 0.13(n=209)
## Avg_parahippocampal_thickavg 2.68 ± 0.22(n=209)
## Avg_paracentral_thickavg 2.49 ± 0.14(n=209)
## Avg_parsopercularis_thickavg 2.65 ± 0.13(n=209)
## Avg_parsorbitalis_thickavg 2.77 ± 0.18(n=209)
## Avg_parstriangularis_thickavg 2.55 ± 0.14(n=209)
## Avg_pericalcarine_thickavg 1.6 ± 0.11(n=209)
## Avg_postcentral_thickavg 2.1 ± 0.11(n=209)
## Avg_posteriorcingulate_thickavg 2.55 ± 0.13(n=209)
## Avg_precentral_thickavg 2.6 ± 0.13(n=209)
## Avg_precuneus_thickavg 2.43 ± 0.12(n=209)
## Avg_rostralanteriorcingulate_thickavg 2.92 ± 0.17(n=209)
## Avg_rostralmiddlefrontal_thickavg 2.48 ± 0.12(n=209)
## Avg_superiorfrontal_thickavg 2.86 ± 0.14(n=209)
## Avg_superiorparietal_thickavg 2.26 ± 0.11(n=209)
## Avg_superiortemporal_thickavg 2.85 ± 0.14(n=209)
## Avg_supramarginal_thickavg 2.6 ± 0.13(n=209)
## Avg_frontalpole_thickavg 2.95 ± 0.24(n=209)
## Avg_temporalpole_thickavg 3.73 ± 0.22(n=209)
## Avg_transversetemporal_thickavg 2.3 ± 0.19(n=209)
## Avg_insula_thickavg 3 ± 0.16(n=209)
## Avg_bankssts_surfavg 944.64 ± 135.77(n=209)
## Avg_caudalanteriorcingulate_surfavg 631.71 ± 121.04(n=209)
## Avg_caudalmiddlefrontal_surfavg 2108.21 ± 324.36(n=209)
## Avg_cuneus_surfavg 1595.19 ± 209.57(n=209)
## Avg_entorhinal_surfavg 417.4 ± 68.72(n=209)
## Avg_fusiform_surfavg 3070.39 ± 380(n=209)
## Avg_inferiorparietal_surfavg 5090.6 ± 691.75(n=209)
## Avg_inferiortemporal_surfavg 3422.45 ± 477.41(n=209)
## Avg_isthmuscingulate_surfavg 985.27 ± 146.75(n=209)
## Avg_lateraloccipital_surfavg 5010.8 ± 641.02(n=209)
## Avg_lateralorbitofrontal_surfavg 2545.06 ± 280.94(n=209)
## Avg_lingual_surfavg 3176.06 ± 405.42(n=209)
## Avg_medialorbitofrontal_surfavg 1889.05 ± 214.34(n=209)
## Avg_middletemporal_surfavg 3389.74 ± 457.63(n=209)
## Avg_parahippocampal_surfavg 652.13 ± 78.2(n=209)
## Avg_paracentral_surfavg 1418.23 ± 164.58(n=209)
## Avg_parsopercularis_surfavg 1474.17 ± 207.84(n=209)
## Avg_parsorbitalis_surfavg 713.9 ± 89.72(n=209)
## Avg_parstriangularis_surfavg 1414.59 ± 212.28(n=209)
## Avg_pericalcarine_surfavg 1544.01 ± 254.77(n=209)
## Avg_postcentral_surfavg 4132.66 ± 475.47(n=209)
## Avg_posteriorcingulate_surfavg 1175.53 ± 170.86(n=209)
## Avg_precentral_surfavg 4904.58 ± 519.49(n=209)
## Avg_precuneus_surfavg 3983.16 ± 522.64(n=209)
## Avg_rostralanteriorcingulate_surfavg 691.98 ± 130.93(n=209)
## Avg_rostralmiddlefrontal_surfavg 5501.88 ± 839.1(n=209)
## Avg_superiorfrontal_surfavg 6706.45 ± 847.65(n=209)
## Avg_superiorparietal_surfavg 5380.33 ± 666.04(n=209)
## Avg_superiortemporal_surfavg 3799.53 ± 458.77(n=209)
## Avg_supramarginal_surfavg 3868.51 ± 590.94(n=209)
## Avg_frontalpole_surfavg 272.64 ± 29.73(n=209)
## Avg_temporalpole_surfavg 470.19 ± 50.68(n=209)
## Avg_transversetemporal_surfavg 405.54 ± 56.68(n=209)
## Avg_insula_surfavg 2427.76 ± 267.23(n=209)
## ICV 1479141.93 ± 178127.94(n=209)
## Avg_LatVent 7859.21 ± 4663.11(n=209)
## Avg_thal 7848.64 ± 875.31(n=209)
## Avg_caud 3460 ± 458.02(n=209)
## Avg_put 4990.17 ± 616.09(n=209)
## Avg_pal 1983.01 ± 233.1(n=209)
## Avg_hippo 4188.57 ± 399.62(n=209)
## Avg_amyg 1720.7 ± 204.65(n=209)
## Avg_accumb 469.26 ± 84.9(n=209)
## Schizophrenia..CNUH. p.value
## Avg_bankssts_thickavg 2.52 ± 0.12(n=196) 0.0000
## Avg_caudalanteriorcingulate_thickavg 2.65 ± 0.19(n=196) 0.0016
## Avg_caudalmiddlefrontal_thickavg 2.54 ± 0.12(n=196) 0.0000
## Avg_cuneus_thickavg 1.86 ± 0.11(n=196) 0.0000
## Avg_entorhinal_thickavg 3.4 ± 0.25(n=196) 0.0106
## Avg_fusiform_thickavg 2.71 ± 0.12(n=196) 0.0000
## Avg_inferiorparietal_thickavg 2.47 ± 0.11(n=196) 0.0000
## Avg_inferiortemporal_thickavg 2.78 ± 0.12(n=196) 0.0000
## Avg_isthmuscingulate_thickavg 2.27 ± 0.16(n=196) 0.0000
## Avg_lateraloccipital_thickavg 2.22 ± 0.1(n=196) 0.0000
## Avg_lateralorbitofrontal_thickavg 2.6 ± 0.13(n=196) 0.0000
## Avg_lingual_thickavg 1.98 ± 0.1(n=196) 0.0011
## Avg_medialorbitofrontal_thickavg 2.47 ± 0.14(n=196) 0.0000
## Avg_middletemporal_thickavg 2.86 ± 0.13(n=196) 0.0000
## Avg_parahippocampal_thickavg 2.6 ± 0.23(n=196) 0.0000
## Avg_paracentral_thickavg 2.44 ± 0.12(n=196) 0.0000
## Avg_parsopercularis_thickavg 2.55 ± 0.14(n=196) 0.0000
## Avg_parsorbitalis_thickavg 2.64 ± 0.16(n=196) 0.0000
## Avg_parstriangularis_thickavg 2.45 ± 0.12(n=196) 0.0000
## Avg_pericalcarine_thickavg 1.58 ± 0.1(n=196) 0.0027
## Avg_postcentral_thickavg 2.03 ± 0.09(n=196) 0.0000
## Avg_posteriorcingulate_thickavg 2.48 ± 0.14(n=196) 0.0000
## Avg_precentral_thickavg 2.52 ± 0.13(n=196) 0.0000
## Avg_precuneus_thickavg 2.36 ± 0.11(n=196) 0.0000
## Avg_rostralanteriorcingulate_thickavg 2.84 ± 0.18(n=196) 0.0000
## Avg_rostralmiddlefrontal_thickavg 2.37 ± 0.1(n=196) 0.0000
## Avg_superiorfrontal_thickavg 2.76 ± 0.12(n=196) 0.0000
## Avg_superiorparietal_thickavg 2.2 ± 0.1(n=196) 0.0000
## Avg_superiortemporal_thickavg 2.74 ± 0.14(n=196) 0.0000
## Avg_supramarginal_thickavg 2.52 ± 0.12(n=196) 0.0000
## Avg_frontalpole_thickavg 2.79 ± 0.23(n=196) 0.0000
## Avg_temporalpole_thickavg 3.62 ± 0.24(n=196) 0.0000
## Avg_transversetemporal_thickavg 2.28 ± 0.17(n=196) 0.0002
## Avg_insula_thickavg 2.94 ± 0.14(n=196) 0.0000
## Avg_bankssts_surfavg 948.13 ± 137.1(n=196) 0.0406
## Avg_caudalanteriorcingulate_surfavg 597.62 ± 107.91(n=196) 0.0000
## Avg_caudalmiddlefrontal_surfavg 2105.85 ± 317.56(n=196) 0.0163
## Avg_cuneus_surfavg 1637.13 ± 211.18(n=196) 0.0784
## Avg_entorhinal_surfavg 417.58 ± 66.72(n=196) 0.9069
## Avg_fusiform_surfavg 3001.54 ± 360.22(n=196) 0.0000
## Avg_inferiorparietal_surfavg 5120.73 ± 659.51(n=196) 0.0314
## Avg_inferiortemporal_surfavg 3437.63 ± 457.53(n=196) 0.0812
## Avg_isthmuscingulate_surfavg 1011.05 ± 139.99(n=196) 0.0860
## Avg_lateraloccipital_surfavg 5050.99 ± 633.71(n=196) 0.0716
## Avg_lateralorbitofrontal_surfavg 2600.96 ± 291.54(n=196) 0.0703
## Avg_lingual_surfavg 3170.42 ± 396.38(n=196) 0.0821
## Avg_medialorbitofrontal_surfavg 1902.11 ± 215.58(n=196) 0.2893
## Avg_middletemporal_surfavg 3416.84 ± 442.17(n=196) 0.0238
## Avg_parahippocampal_surfavg 642.94 ± 75.34(n=196) 0.0040
## Avg_paracentral_surfavg 1443.73 ± 171.59(n=196) 0.2557
## Avg_parsopercularis_surfavg 1452.61 ± 223.15(n=196) 0.0054
## Avg_parsorbitalis_surfavg 733.48 ± 96.98(n=196) 0.0528
## Avg_parstriangularis_surfavg 1425.49 ± 224.21(n=196) 0.0619
## Avg_pericalcarine_surfavg 1553.84 ± 242.96(n=196) 0.5293
## Avg_postcentral_surfavg 4163 ± 469.17(n=196) 0.1896
## Avg_posteriorcingulate_surfavg 1167.49 ± 169.67(n=196) 0.0107
## Avg_precentral_surfavg 4967.66 ± 519.56(n=196) 0.1567
## Avg_precuneus_surfavg 4018.58 ± 496.78(n=196) 0.0542
## Avg_rostralanteriorcingulate_surfavg 678.14 ± 129.5(n=196) 0.0079
## Avg_rostralmiddlefrontal_surfavg 5744.03 ± 864.2(n=196) 0.0054
## Avg_superiorfrontal_surfavg 6819.75 ± 835.58(n=196) 0.0648
## Avg_superiorparietal_surfavg 5442.31 ± 631.79(n=196) 0.0710
## Avg_superiortemporal_surfavg 3800.74 ± 431.61(n=196) 0.0177
## Avg_supramarginal_surfavg 3906.44 ± 564.67(n=196) 0.1533
## Avg_frontalpole_surfavg 276.63 ± 31.96(n=196) 0.1308
## Avg_temporalpole_surfavg 479.5 ± 54.55(n=196) 0.1792
## Avg_transversetemporal_surfavg 394.04 ± 54.84(n=196) 0.0001
## Avg_insula_surfavg 2416.67 ± 262.66(n=196) 0.0750
## ICV 1508777.92 ± 155869.81(n=196) 0.0526
## Avg_LatVent 9855 ± 5140.15(n=196) 0.0000
## Avg_thal 7204.15 ± 879.25(n=196) 0.0000
## Avg_caud 3391.64 ± 425.19(n=196) 0.0000
## Avg_put 4971.56 ± 637.67(n=196) 0.0001
## Avg_pal 2096.43 ± 234.23(n=196) 0.0000
## Avg_hippo 3945.11 ± 435.22(n=196) 0.0000
## Avg_amyg 1664.91 ± 230.46(n=196) 0.0000
## Avg_accumb 455.59 ± 84.96(n=196) 0.0001
## Group.1.vs.Group.2 Group.1.vs.Group.3
## Avg_bankssts_thickavg 0.1927 0
## Avg_caudalanteriorcingulate_thickavg 0.4142 0.0012
## Avg_caudalmiddlefrontal_thickavg 0.0861 0
## Avg_cuneus_thickavg 0.1531 0
## Avg_entorhinal_thickavg 1 0.0082
## Avg_fusiform_thickavg 0.0721 0
## Avg_inferiorparietal_thickavg 0.0456 0
## Avg_inferiortemporal_thickavg 0.3127 0
## Avg_isthmuscingulate_thickavg 0.2106 0
## Avg_lateraloccipital_thickavg 0.3895 0
## Avg_lateralorbitofrontal_thickavg 0.0178 0
## Avg_lingual_thickavg 0.0761 0.0018
## Avg_medialorbitofrontal_thickavg 0.036 0
## Avg_middletemporal_thickavg 0.08 0
## Avg_parahippocampal_thickavg 1 0
## Avg_paracentral_thickavg 0.2343 0
## Avg_parsopercularis_thickavg 0.0115 0
## Avg_parsorbitalis_thickavg 0.0928 0
## Avg_parstriangularis_thickavg 0.0363 0
## Avg_pericalcarine_thickavg 0.2684 0.0025
## Avg_postcentral_thickavg 0.255 0
## Avg_posteriorcingulate_thickavg 0.0478 0
## Avg_precentral_thickavg 0.1142 0
## Avg_precuneus_thickavg 0.0431 0
## Avg_rostralanteriorcingulate_thickavg 0.149 0
## Avg_rostralmiddlefrontal_thickavg 0.1316 0
## Avg_superiorfrontal_thickavg 0.0561 0
## Avg_superiorparietal_thickavg 0.2885 0
## Avg_superiortemporal_thickavg 0.0112 0
## Avg_supramarginal_thickavg 0.0621 0
## Avg_frontalpole_thickavg 0.5615 0
## Avg_temporalpole_thickavg 0.8051 0
## Avg_transversetemporal_thickavg 0.0783 2e-04
## Avg_insula_thickavg 0.0017 0
## Avg_bankssts_surfavg 0.0851 0.208
## Avg_caudalanteriorcingulate_surfavg 0.3857 0
## Avg_caudalmiddlefrontal_surfavg 0.0696 0.0629
## Avg_cuneus_surfavg 0.1763 1
## Avg_entorhinal_surfavg 1 1
## Avg_fusiform_surfavg 0.1083 0
## Avg_inferiorparietal_surfavg 0.0544 0.2365
## Avg_inferiortemporal_surfavg 0.1402 0.3737
## Avg_isthmuscingulate_surfavg 0.1196 1
## Avg_lateraloccipital_surfavg 0.0908 0.5583
## Avg_lateralorbitofrontal_surfavg 0.1367 1
## Avg_lingual_surfavg 0.2587 0.194
## Avg_medialorbitofrontal_surfavg 0.3689 1
## Avg_middletemporal_surfavg 0.0348 0.2672
## Avg_parahippocampal_surfavg 0.2468 0.0041
## Avg_paracentral_surfavg 0.7685 1
## Avg_parsopercularis_surfavg 0.195 0.0067
## Avg_parsorbitalis_surfavg 0.13 1
## Avg_parstriangularis_surfavg 0.0898 0.4171
## Avg_pericalcarine_surfavg 0.8308 1
## Avg_postcentral_surfavg 0.2316 1
## Avg_posteriorcingulate_surfavg 0.1078 0.024
## Avg_precentral_surfavg 0.1672 1
## Avg_precuneus_surfavg 0.066 0.5214
## Avg_rostralanteriorcingulate_surfavg 0.3032 0.0083
## Avg_rostralmiddlefrontal_surfavg 0.0347 0.6262
## Avg_superiorfrontal_surfavg 0.0584 1
## Avg_superiorparietal_surfavg 0.0709 0.9356
## Avg_superiortemporal_surfavg 0.0628 0.079
## Avg_supramarginal_surfavg 0.1851 0.9545
## Avg_frontalpole_surfavg 0.1386 1
## Avg_temporalpole_surfavg 1 0.5736
## Avg_transversetemporal_surfavg 0.2938 1e-04
## Avg_insula_surfavg 0.3751 0.1252
## ICV 0.0543 1
## Avg_LatVent 1 0
## Avg_thal 0.0039 0
## Avg_caud 0.0312 1e-04
## Avg_put 0.0041 0.0014
## Avg_pal 0.1387 3e-04
## Avg_hippo 0.0383 0
## Avg_amyg 0.0632 0
## Avg_accumb 0.0732 1e-04
## Group.2.vs.Group.3 age.sex.edu.adj.p
## Avg_bankssts_thickavg 0 0
## Avg_caudalanteriorcingulate_thickavg 0.2314 0.0394
## Avg_caudalmiddlefrontal_thickavg 0 0
## Avg_cuneus_thickavg 0.0326 0.0027
## Avg_entorhinal_thickavg 0.1047 0.0086
## Avg_fusiform_thickavg 0 0
## Avg_inferiorparietal_thickavg 0 0
## Avg_inferiortemporal_thickavg 0 0
## Avg_isthmuscingulate_thickavg 0 0
## Avg_lateraloccipital_thickavg 0 0
## Avg_lateralorbitofrontal_thickavg 0 0
## Avg_lingual_thickavg 0.8757 0.1268
## Avg_medialorbitofrontal_thickavg 0 0
## Avg_middletemporal_thickavg 0 0
## Avg_parahippocampal_thickavg 0.0015 0.0012
## Avg_paracentral_thickavg 1e-04 0
## Avg_parsopercularis_thickavg 0 0
## Avg_parsorbitalis_thickavg 0 0
## Avg_parstriangularis_thickavg 0 0
## Avg_pericalcarine_thickavg 0.4702 0.0388
## Avg_postcentral_thickavg 0 0
## Avg_posteriorcingulate_thickavg 0 0
## Avg_precentral_thickavg 0 0
## Avg_precuneus_thickavg 0 0
## Avg_rostralanteriorcingulate_thickavg 0 0
## Avg_rostralmiddlefrontal_thickavg 0 0
## Avg_superiorfrontal_thickavg 0 0
## Avg_superiorparietal_thickavg 0 0
## Avg_superiortemporal_thickavg 0 0
## Avg_supramarginal_thickavg 0 0
## Avg_frontalpole_thickavg 0 0
## Avg_temporalpole_thickavg 0 0
## Avg_transversetemporal_thickavg 0.3634 0.1065
## Avg_insula_thickavg 0 0
## Avg_bankssts_surfavg 1 0.7081
## Avg_caudalanteriorcingulate_surfavg 0.0063 0
## Avg_caudalmiddlefrontal_surfavg 1 0.5272
## Avg_cuneus_surfavg 0.1043 0.3277
## Avg_entorhinal_surfavg 1 0.6574
## Avg_fusiform_surfavg 0.1348 9e-04
## Avg_inferiorparietal_surfavg 1 0.7843
## Avg_inferiortemporal_surfavg 1 0.461
## Avg_isthmuscingulate_surfavg 0.1799 0.441
## Avg_lateraloccipital_surfavg 1 0.7346
## Avg_lateralorbitofrontal_surfavg 0.1056 0.1934
## Avg_lingual_surfavg 1 0.5824
## Avg_medialorbitofrontal_surfavg 1 0.4816
## Avg_middletemporal_surfavg 1 0.6729
## Avg_parahippocampal_surfavg 0.6169 0.0469
## Avg_paracentral_surfavg 0.3043 0.6717
## Avg_parsopercularis_surfavg 0.8696 0.1417
## Avg_parsorbitalis_surfavg 0.0684 0.2962
## Avg_parstriangularis_surfavg 1 0.8697
## Avg_pericalcarine_surfavg 1 0.9883
## Avg_postcentral_surfavg 1 0.5574
## Avg_posteriorcingulate_surfavg 1 0.2626
## Avg_precentral_surfavg 0.5791 0.7428
## Avg_precuneus_surfavg 1 0.8551
## Avg_rostralanteriorcingulate_surfavg 0.725 0.0086
## Avg_rostralmiddlefrontal_surfavg 0.0055 0.013
## Avg_superiorfrontal_surfavg 0.422 0.5975
## Avg_superiorparietal_surfavg 0.941 0.7149
## Avg_superiortemporal_surfavg 1 0.337
## Avg_supramarginal_surfavg 1 0.5834
## Avg_frontalpole_surfavg 0.4631 0.756
## Avg_temporalpole_surfavg 0.197 0.5612
## Avg_transversetemporal_surfavg 0.0946 0.0021
## Avg_insula_surfavg 1 0.0757
## ICV 0.1968 0.7287
## Avg_LatVent 0 2e-04
## Avg_thal 0 0
## Avg_caud 0.3826 0.0188
## Avg_put 1 0.1725
## Avg_pal 0 0
## Avg_hippo 0 0
## Avg_amyg 0.0164 0
## Avg_accumb 0.2939 0.0614
## bonf.adj.p.value
## Avg_bankssts_thickavg 0.0000
## Avg_caudalanteriorcingulate_thickavg 0.1232
## Avg_caudalmiddlefrontal_thickavg 0.0000
## Avg_cuneus_thickavg 0.0000
## Avg_entorhinal_thickavg 0.8162
## Avg_fusiform_thickavg 0.0000
## Avg_inferiorparietal_thickavg 0.0000
## Avg_inferiortemporal_thickavg 0.0000
## Avg_isthmuscingulate_thickavg 0.0000
## Avg_lateraloccipital_thickavg 0.0000
## Avg_lateralorbitofrontal_thickavg 0.0000
## Avg_lingual_thickavg 0.0847
## Avg_medialorbitofrontal_thickavg 0.0000
## Avg_middletemporal_thickavg 0.0000
## Avg_parahippocampal_thickavg 0.0000
## Avg_paracentral_thickavg 0.0000
## Avg_parsopercularis_thickavg 0.0000
## Avg_parsorbitalis_thickavg 0.0000
## Avg_parstriangularis_thickavg 0.0000
## Avg_pericalcarine_thickavg 0.2079
## Avg_postcentral_thickavg 0.0000
## Avg_posteriorcingulate_thickavg 0.0000
## Avg_precentral_thickavg 0.0000
## Avg_precuneus_thickavg 0.0000
## Avg_rostralanteriorcingulate_thickavg 0.0000
## Avg_rostralmiddlefrontal_thickavg 0.0000
## Avg_superiorfrontal_thickavg 0.0000
## Avg_superiorparietal_thickavg 0.0000
## Avg_superiortemporal_thickavg 0.0000
## Avg_supramarginal_thickavg 0.0000
## Avg_frontalpole_thickavg 0.0000
## Avg_temporalpole_thickavg 0.0000
## Avg_transversetemporal_thickavg 0.0154
## Avg_insula_thickavg 0.0000
## Avg_bankssts_surfavg 1.0000
## Avg_caudalanteriorcingulate_surfavg 0.0000
## Avg_caudalmiddlefrontal_surfavg 1.0000
## Avg_cuneus_surfavg 1.0000
## Avg_entorhinal_surfavg 1.0000
## Avg_fusiform_surfavg 0.0000
## Avg_inferiorparietal_surfavg 1.0000
## Avg_inferiortemporal_surfavg 1.0000
## Avg_isthmuscingulate_surfavg 1.0000
## Avg_lateraloccipital_surfavg 1.0000
## Avg_lateralorbitofrontal_surfavg 1.0000
## Avg_lingual_surfavg 1.0000
## Avg_medialorbitofrontal_surfavg 1.0000
## Avg_middletemporal_surfavg 1.0000
## Avg_parahippocampal_surfavg 0.3080
## Avg_paracentral_surfavg 1.0000
## Avg_parsopercularis_surfavg 0.4158
## Avg_parsorbitalis_surfavg 1.0000
## Avg_parstriangularis_surfavg 1.0000
## Avg_pericalcarine_surfavg 1.0000
## Avg_postcentral_surfavg 1.0000
## Avg_posteriorcingulate_surfavg 0.8239
## Avg_precentral_surfavg 1.0000
## Avg_precuneus_surfavg 1.0000
## Avg_rostralanteriorcingulate_surfavg 0.6083
## Avg_rostralmiddlefrontal_surfavg 0.4158
## Avg_superiorfrontal_surfavg 1.0000
## Avg_superiorparietal_surfavg 1.0000
## Avg_superiortemporal_surfavg 1.0000
## Avg_supramarginal_surfavg 1.0000
## Avg_frontalpole_surfavg 1.0000
## Avg_temporalpole_surfavg 1.0000
## Avg_transversetemporal_surfavg 0.0077
## Avg_insula_surfavg 1.0000
## ICV 1.0000
## Avg_LatVent 0.0000
## Avg_thal 0.0000
## Avg_caud 0.0000
## Avg_put 0.0077
## Avg_pal 0.0000
## Avg_hippo 0.0000
## Avg_amyg 0.0000
## Avg_accumb 0.0077
## bonf.age.sex.edu.adj.p.value
## Avg_bankssts_thickavg 0.0000
## Avg_caudalanteriorcingulate_thickavg 1.0000
## Avg_caudalmiddlefrontal_thickavg 0.0000
## Avg_cuneus_thickavg 0.2079
## Avg_entorhinal_thickavg 0.6622
## Avg_fusiform_thickavg 0.0000
## Avg_inferiorparietal_thickavg 0.0000
## Avg_inferiortemporal_thickavg 0.0000
## Avg_isthmuscingulate_thickavg 0.0000
## Avg_lateraloccipital_thickavg 0.0000
## Avg_lateralorbitofrontal_thickavg 0.0000
## Avg_lingual_thickavg 1.0000
## Avg_medialorbitofrontal_thickavg 0.0000
## Avg_middletemporal_thickavg 0.0000
## Avg_parahippocampal_thickavg 0.0924
## Avg_paracentral_thickavg 0.0000
## Avg_parsopercularis_thickavg 0.0000
## Avg_parsorbitalis_thickavg 0.0000
## Avg_parstriangularis_thickavg 0.0000
## Avg_pericalcarine_thickavg 1.0000
## Avg_postcentral_thickavg 0.0000
## Avg_posteriorcingulate_thickavg 0.0000
## Avg_precentral_thickavg 0.0000
## Avg_precuneus_thickavg 0.0000
## Avg_rostralanteriorcingulate_thickavg 0.0000
## Avg_rostralmiddlefrontal_thickavg 0.0000
## Avg_superiorfrontal_thickavg 0.0000
## Avg_superiorparietal_thickavg 0.0000
## Avg_superiortemporal_thickavg 0.0000
## Avg_supramarginal_thickavg 0.0000
## Avg_frontalpole_thickavg 0.0000
## Avg_temporalpole_thickavg 0.0000
## Avg_transversetemporal_thickavg 1.0000
## Avg_insula_thickavg 0.0000
## Avg_bankssts_surfavg 1.0000
## Avg_caudalanteriorcingulate_surfavg 0.0000
## Avg_caudalmiddlefrontal_surfavg 1.0000
## Avg_cuneus_surfavg 1.0000
## Avg_entorhinal_surfavg 1.0000
## Avg_fusiform_surfavg 0.0693
## Avg_inferiorparietal_surfavg 1.0000
## Avg_inferiortemporal_surfavg 1.0000
## Avg_isthmuscingulate_surfavg 1.0000
## Avg_lateraloccipital_surfavg 1.0000
## Avg_lateralorbitofrontal_surfavg 1.0000
## Avg_lingual_surfavg 1.0000
## Avg_medialorbitofrontal_surfavg 1.0000
## Avg_middletemporal_surfavg 1.0000
## Avg_parahippocampal_surfavg 1.0000
## Avg_paracentral_surfavg 1.0000
## Avg_parsopercularis_surfavg 1.0000
## Avg_parsorbitalis_surfavg 1.0000
## Avg_parstriangularis_surfavg 1.0000
## Avg_pericalcarine_surfavg 1.0000
## Avg_postcentral_surfavg 1.0000
## Avg_posteriorcingulate_surfavg 1.0000
## Avg_precentral_surfavg 1.0000
## Avg_precuneus_surfavg 1.0000
## Avg_rostralanteriorcingulate_surfavg 0.6622
## Avg_rostralmiddlefrontal_surfavg 1.0000
## Avg_superiorfrontal_surfavg 1.0000
## Avg_superiorparietal_surfavg 1.0000
## Avg_superiortemporal_surfavg 1.0000
## Avg_supramarginal_surfavg 1.0000
## Avg_frontalpole_surfavg 1.0000
## Avg_temporalpole_surfavg 1.0000
## Avg_transversetemporal_surfavg 0.1617
## Avg_insula_surfavg 1.0000
## ICV 1.0000
## Avg_LatVent 0.0154
## Avg_thal 0.0000
## Avg_caud 1.0000
## Avg_put 1.0000
## Avg_pal 0.0000
## Avg_hippo 0.0000
## Avg_amyg 0.0000
## Avg_accumb 1.0000# add sex in training
Train_HC_data$Sex <- as.numeric(Train_HC_data$Sex)
Test_HC_data$Sex <- as.numeric(Test_HC_data$Sex)
Test_Patient_dat$Sex <- as.numeric(Test_Patient_dat$Sex)
Test_schizophrenia_dat$Sex <- as.numeric(Test_schizophrenia_dat$Sex)
Test_TRS_dat$Sex <- as.numeric(Test_TRS_dat$Sex)
Test_FESSD_dat$Sex <- as.numeric(Test_FESSD_dat$Sex)Ridge (Nested CV)
library(caret)## 필요한 패키지를 로딩중입니다: latticelibrary(parallel)
Nfolds <- 10
f <- as.formula(paste("Age~Sex+",paste(colnames(Train_HC_data)[5:81],collapse = "+"),sep=""))
fitcontrol <- trainControl(method = "repeatedcv", number=Nfolds, search = "random")
ridge.fun <- function(N,i){
Nfolds <- N
set.seed(i)
outer_folds <- createFolds(Train_HC_data$Age,k=Nfolds)
cv.glmnet <- lapply(outer_folds,function(index) {
tr <- Train_HC_data[-index,]
ts <- Train_HC_data[index,]
fit <- train(f,data=tr,method="ridge",metric="MAE",
trControl=fitcontrol)
c(postResample(predict(fit,ts),ts$Age),lambda=fit$bestTune$lambda)
})
}
Ntimes <- 10
RES <- NULL
for(i in 1:Ntimes){
res <- mclapply(10,i,FUN=ridge.fun,mc.cores = 6)
RES <- rbind(RES,t(data.frame(res)))
}
library(glmnet)## 필요한 패키지를 로딩중입니다: Matrix## Loaded glmnet 4.1-6ridge.fit <- glmnet(as.matrix(Train_HC_data[,c(3,5:81)]),Train_HC_data[,2],
alpha=0,lambda = RES[,4][which.min(RES[,3])])
# Train Control
ctrl_pred_age <- predict(ridge.fit, newx = as.matrix(Train_HC_data[,c(3,5:81)]))
train_error_ctrl_rid <- ctrl_pred_age-Train_HC_data$Age
par(mfrow=c(1,1))
plot(ctrl_pred_age,Train_HC_data$Age,
xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (Train control data)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",
legend=c(paste("r=",round(cor(Train_HC_data$Age,ctrl_pred_age),4)),
paste("R2=",round(cor(Train_HC_data$Age,ctrl_pred_age)^2,4)),
paste("MAE=",round(sum(abs(train_error_ctrl_rid)/length(train_error_ctrl_rid)),4)),
paste("VEcv=",round(1-sum(train_error_ctrl_rid^2/(var(Train_HC_data$Age)*(length(Train_HC_data$Age)-1))),4)*100,"%")) ,bty="n")
# Test Control
ctrl_pred_age <- predict(ridge.fit, newx = as.matrix(Test_HC_data[,c(3,5:81)]))
test_error_ctrl_rid <- ctrl_pred_age-Test_HC_data$Age
plot(ctrl_pred_age,Test_HC_data$Age,
xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (Test control data)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",
legend=c(paste("r=",round(cor(Test_HC_data$Age,ctrl_pred_age),4)),
paste("R2=",round(cor(Test_HC_data$Age,ctrl_pred_age)^2,4)),
paste("MAE=",round(sum(abs(test_error_ctrl_rid)/length(test_error_ctrl_rid)),4)),
paste("VEcv=",round(1-sum(test_error_ctrl_rid^2/(var(Test_HC_data$Age)*(length(Test_HC_data$Age)-1))),4)*100,"%")) ,bty="n")
# Patient
pat_pred_age <- predict(ridge.fit, newx = as.matrix(Test_Patient_dat[,c(3,5:81)]))
test_error_pat_rid <- pat_pred_age-Test_Patient_dat$Age
plot(pat_pred_age,Test_Patient_dat$Age,xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (Patient)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",legend=c(paste("r=",round(cor(Test_Patient_dat$Age,pat_pred_age),4)),
paste("R2=",round(cor(Test_Patient_dat$Age,pat_pred_age)^2,4)),
paste("MAE=",round(sum(abs(test_error_pat_rid)/length(test_error_pat_rid)),4)),
paste("VEcv=",round(1-sum(test_error_pat_rid^2/(var(Test_Patient_dat$Age)*(length(Test_Patient_dat$Age)-1))),4)*100,"%")),bty="n")
# Patient - Schizo
sch_pred_age <- predict(ridge.fit, newx = as.matrix(Test_schizophrenia_dat[,c(3,5:81)]))
test_error_sch_rid <- sch_pred_age-Test_schizophrenia_dat$Age
plot(sch_pred_age,Test_schizophrenia_dat$Age,xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (Schizophrenia)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",legend=c(paste("r=",round(cor(Test_schizophrenia_dat$Age,sch_pred_age),4)),
paste("R2=",round(cor(Test_schizophrenia_dat$Age,sch_pred_age)^2,4)),
paste("MAE=",round(sum(abs(test_error_sch_rid)/length(test_error_sch_rid)),4)),
paste("VEcv=",round(1-sum(test_error_sch_rid^2/(var(Test_schizophrenia_dat$Age)*(length(Test_schizophrenia_dat$Age)-1))),4)*100,"%")),bty="n")
# Patient - TRS
TRS_pred_age <- predict(ridge.fit, newx = as.matrix(Test_TRS_dat[,c(3,5:81)]))
test_error_trs_rid <- TRS_pred_age-Test_TRS_dat$Age
plot(TRS_pred_age,Test_TRS_dat$Age,xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (TRS)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",legend=c(paste("r=",round(cor(Test_TRS_dat$Age,TRS_pred_age),4)),
paste("R2=",round(cor(Test_TRS_dat$Age,TRS_pred_age)^2,4)),
paste("MAE=",round(sum(abs(test_error_trs_rid)/length(test_error_trs_rid)),4)),
paste("VEcv=",round(1-sum(test_error_trs_rid^2/(var(Test_TRS_dat$Age)*(length(Test_TRS_dat$Age)-1))),4)*100,"%")),bty="n")
# Patient - FESSD
SSD_pred_age <- predict(ridge.fit, newx = as.matrix(Test_FESSD_dat[,c(3,5:81)]))
test_error_ssd_rid <- SSD_pred_age-Test_FESSD_dat$Age
plot(SSD_pred_age,Test_FESSD_dat$Age,xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (FE-SSD)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",legend=c(paste("r=",round(cor(Test_FESSD_dat$Age,SSD_pred_age),4)),
paste("R2=",round(cor(Test_FESSD_dat$Age,SSD_pred_age)^2,4)),
paste("MAE=",round(sum(abs(test_error_ssd_rid)/length(test_error_ssd_rid)),4)),
paste("VEcv=",round(1-sum(test_error_ssd_rid^2/(var(Test_FESSD_dat$Age)*(length(Test_FESSD_dat$Age)-1))),4)*100,"%")),bty="n")
Brain PAD correction-plot
#### Brain PAD correction-HC
pred_age <- predict(ridge.fit,newx=as.matrix(Test_HC_data[,c(3,5:81)]))
test_error_rid <- pred_age-Test_HC_data$Age
PAD_rid <- test_error_rid
cor.test(PAD_rid,Test_HC_data$Age)$p.value## [1] 2.735224e-17png("1_rigde_corr.png", width = 800, height = 600)
par(mfrow=c(2,2), oma=c(1,3,1,1))
plot(Test_HC_data$Age,PAD_rid,xlim=c(10,80),ylim=c(-40,40),pch=16,
xlab="Chronological age (Test HC)",ylab="Brain PAD uncorrected")
abline(h=0,lty=2,col="red")
abline(lm(PAD_rid~Test_HC_data$Age))
legend("topright",legend=paste("r=",round(cor(PAD_rid,Test_HC_data$Age),4),"(p<0.05)",sep=""))
mtext("A", side=3, line=-3, cex=2, adj = -0.0, outer=TRUE)
Model_linear <- lm(PAD_rid~Test_HC_data$Age+Test_HC_data$Sex)
NullModel <- lm(PAD_rid~1)
anova(NullModel,Model_linear,test="Chisq") # Reject H0 (Full model is better)## Analysis of Variance Table
##
## Model 1: PAD_rid ~ 1
## Model 2: PAD_rid ~ Test_HC_data$Age + Test_HC_data$Sex
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 208 16868
## 2 206 11881 2 4987.8 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Test_HC_data$Age2 <- Test_HC_data$Age^2
Model_quad <- lm(PAD_rid~Test_HC_data$Sex+Test_HC_data$Age+Test_HC_data$Age2)
anova(NullModel,Model_quad,test="Chisq") # Reject H0 (Full model is better)## Analysis of Variance Table
##
## Model 1: PAD_rid ~ 1
## Model 2: PAD_rid ~ Test_HC_data$Sex + Test_HC_data$Age + Test_HC_data$Age2
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 208 16868
## 2 205 11818 3 5050.3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(Model_linear,Model_quad,test="Chisq") # Not Reject H0 (Reduced model is better)## Analysis of Variance Table
##
## Model 1: PAD_rid ~ Test_HC_data$Age + Test_HC_data$Sex
## Model 2: PAD_rid ~ Test_HC_data$Sex + Test_HC_data$Age + Test_HC_data$Age2
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 206 11881
## 2 205 11818 1 62.413 0.2981cPAD_rid <- lm(PAD_rid~+Test_HC_data$Sex+Test_HC_data$Age)$residuals
plot(Test_HC_data$Age,cPAD_rid,xlim=c(10,75),ylim=c(-40,40),pch=16,
xlab="Chronological age (Test HC)",ylab="Brain PAD corrected for age")
abline(lm(cPAD_rid~Test_HC_data$Age))
abline(h=0,lty=2,col="red")
mtext("B", side=3, line=-3, cex=2, adj = 0.5, outer=TRUE)
#### Brain PAD correction-Patient
pred_age <- predict(ridge.fit,newx=as.matrix(Test_Patient_dat[,c(3,5:81)]))
test_error_rid <- pred_age-Test_Patient_dat$Age
PAD_rid <- test_error_rid
plot(Test_Patient_dat$Age,PAD_rid,xlim=c(10,80),ylim=c(-40,40),pch=16,
xlab="Chronological age (Test Patient)",ylab="Brain PAD uncorrected")
abline(h=0,lty=2,col="red")
abline(lm(PAD_rid~Test_Patient_dat$Age))
legend("topright",legend=paste("r=",round(cor(PAD_rid,Test_Patient_dat$Age),4),"(p<0.05)",sep=""))
cor.test(PAD_rid,Test_Patient_dat$Age)$p.value## [1] 3.574547e-10mtext("C", side=3, line=-27, cex=2, adj = -0.0, outer=TRUE)
Model_linear <- lm(PAD_rid~Test_Patient_dat$Age+Test_Patient_dat$Sex)
NullModel <- lm(PAD_rid~1)
anova(NullModel,Model_linear,test="Chisq") # Reject H0 (Full model is better)## Analysis of Variance Table
##
## Model 1: PAD_rid ~ 1
## Model 2: PAD_rid ~ Test_Patient_dat$Age + Test_Patient_dat$Sex
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 232 24702
## 2 230 20756 2 3946 3.2e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Test_Patient_dat$Age2 <- Test_Patient_dat$Age^2
Model_quad <- lm(PAD_rid~Test_Patient_dat$Sex+Test_Patient_dat$Age+Test_Patient_dat$Age2)
anova(NullModel,Model_quad,test="Chisq") # Reject H0 (Full model is better)## Analysis of Variance Table
##
## Model 1: PAD_rid ~ 1
## Model 2: PAD_rid ~ Test_Patient_dat$Sex + Test_Patient_dat$Age + Test_Patient_dat$Age2
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 232 24702
## 2 229 20335 3 4367.5 1.191e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(Model_linear,Model_quad,test="Chisq") # Reject H0 (Full model is better)## Analysis of Variance Table
##
## Model 1: PAD_rid ~ Test_Patient_dat$Age + Test_Patient_dat$Sex
## Model 2: PAD_rid ~ Test_Patient_dat$Sex + Test_Patient_dat$Age + Test_Patient_dat$Age2
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 230 20756
## 2 229 20335 1 421.53 0.02935 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cPAD_rid <- lm(PAD_rid~+Test_Patient_dat$Sex+Test_Patient_dat$Age+Test_Patient_dat$Age2)$residuals
plot(Test_Patient_dat$Age,cPAD_rid,xlim=c(10,75),ylim=c(-40,40),pch=16,
xlab="Chronological age (Test Patient)",ylab="Brain PAD corrected for age")
abline(lm(cPAD_rid~Test_Patient_dat$Age))
abline(h=0,lty=2,col="red")
mtext("D", side=3, line=-27, cex=2, adj = 0.5, outer=TRUE)
dev.off()## png
## 2Brain PAD correction -table
test_dat <- rbind(Test_HC_data,Test_Patient_dat)
test_dat$Dx <- c(rep(0,nrow(Test_HC_data)),rep(1,nrow(Test_Patient_dat)))
test_dat$Dx <- factor(test_dat$Dx)
pred_age <- predict(ridge.fit,newx=as.matrix(test_dat[,c(3,5:81)]))
test_error_rid <- pred_age-test_dat$Age
PAD_rid <- test_error_rid
mean(PAD_rid) ; sd(PAD_rid) ; median(PAD_rid)## [1] 1.513638## [1] 10.23045## [1] 2.005724Model_linear <- lm(PAD_rid~test_dat$Age+test_dat$Sex)
NullModel <- lm(PAD_rid~1)
anova(NullModel,Model_linear,test="Chisq") # Reject H0 (Full model is better)## Analysis of Variance Table
##
## Model 1: PAD_rid ~ 1
## Model 2: PAD_rid ~ test_dat$Age + test_dat$Sex
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 441 46156
## 2 439 36892 2 9264.3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1test_dat$Age2 <- test_dat$Age^2
Model_quad <- lm(PAD_rid~test_dat$Sex+test_dat$Age+test_dat$Age2)
anova(NullModel,Model_quad,test="Chisq") # Reject H0 (Full model is better)## Analysis of Variance Table
##
## Model 1: PAD_rid ~ 1
## Model 2: PAD_rid ~ test_dat$Sex + test_dat$Age + test_dat$Age2
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 441 46156
## 2 438 36157 3 9999.4 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(Model_linear,Model_quad,test="Chisq") # Reject H0 ## Analysis of Variance Table
##
## Model 1: PAD_rid ~ test_dat$Age + test_dat$Sex
## Model 2: PAD_rid ~ test_dat$Sex + test_dat$Age + test_dat$Age2
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 439 36892
## 2 438 36157 1 735.1 0.002844 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cPAD_rid <- lm(PAD_rid~+test_dat$Sex+test_dat$Age+test_dat$Age2)$residuals
mean(cPAD_rid) ; sd(cPAD_rid) ; median(cPAD_rid)## [1] -4.443404e-16## [1] 9.054699## [1] 0.02024823library(plotrix)
dat <- data.frame(SubjID=test_dat$SubjID,Brain_PAD = c(PAD_rid),Status=rep(c("Control","Patient"),c(nrow(Test_HC_data),nrow(Test_Patient_dat))))
dat$status_sch <- dat$status_trs <- dat$status_ssd <- rep(NA,nrow(dat))
dat$status_sch[test_dat$SubjID %in% Test_HC_data$SubjID ] <- "Control"
dat$status_sch[test_dat$SubjID %in% Test_schizophrenia_dat$SubjID ] <- "Schizophrenia"
dat$status_trs[test_dat$SubjID %in% Test_TRS_dat$SubjID ] <- "TRS"
dat$status_trs[test_dat$SubjID %in% Test_HC_data$SubjID ] <- "Control"
dat$status_ssd[test_dat$SubjID %in% Test_FESSD_dat$SubjID ] <- "FE-SSD"
dat$status_ssd[test_dat$SubjID %in% Test_HC_data$SubjID ] <- "Control"
library(dplyr)
palette <- c("FE-SSD"="darkgreen","TRS"="red","Control"="grey","Schizophrenia"="blue","Patient"="orange")
dat %>% dplyr::select(Brain_PAD,Status) %>% group_by(Status) %>% summarise_each(funs(mean,sd,std.error))## Warning: `summarise_each_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `across()` instead.
## ℹ The deprecated feature was likely used in the dplyr package.
## Please report the issue at <https://github.com/tidyverse/dplyr/issues>.## Warning: `funs()` was deprecated in dplyr 0.8.0.
## ℹ Please use a list of either functions or lambdas:
##
## # Simple named list: list(mean = mean, median = median)
##
## # Auto named with `tibble::lst()`: tibble::lst(mean, median)
##
## # Using lambdas list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))## # A tibble: 2 × 4
## Status mean sd std.error
## <chr> <dbl> <dbl> <dbl>
## 1 Control -1.89 9.01 0.623
## 2 Patient 4.56 10.3 0.676dat %>% dplyr::select(Brain_PAD,status_sch) %>% group_by(status_sch) %>% summarise_each(funs(mean,sd,std.error))## # A tibble: 3 × 4
## status_sch mean sd std.error
## <chr> <dbl> <dbl> <dbl>
## 1 Control -1.89 9.01 0.623
## 2 Schizophrenia 4.58 10.5 0.750
## 3 <NA> 4.47 9.42 1.55dat %>% dplyr::select(Brain_PAD,status_trs) %>% group_by(status_trs) %>% summarise_each(funs(mean,sd,std.error))## # A tibble: 3 × 4
## status_trs mean sd std.error
## <chr> <dbl> <dbl> <dbl>
## 1 Control -1.89 9.01 0.623
## 2 TRS 6.97 7.54 1.11
## 3 <NA> 3.97 10.8 0.792dat %>% dplyr::select(Brain_PAD,status_ssd) %>% group_by(status_ssd) %>% summarise_each(funs(mean,sd,std.error))## # A tibble: 3 × 4
## status_ssd mean sd std.error
## <chr> <dbl> <dbl> <dbl>
## 1 Control -1.89 9.01 0.623
## 2 FE-SSD 5.19 10.7 1.26
## 3 <NA> 4.29 10.2 0.801library(ggplot2)
ggplot(dat[!is.na(dat$status_sch),],aes(x = Brain_PAD, fill = status_sch)) + scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))
ggplot(dat[!is.na(dat$status_ssd),],aes(x = Brain_PAD, fill = status_ssd)) + scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))
ggplot(dat[!is.na(dat$status_trs),],aes(x = Brain_PAD, fill = status_trs)) +scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))
summary(lm(Brain_PAD~Status+test_dat$Sex+test_dat$Age+test_dat$Age2,data=dat))##
## Call:
## lm(formula = Brain_PAD ~ Status + test_dat$Sex + test_dat$Age +
## test_dat$Age2, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -41.862 -5.379 0.325 5.900 23.850
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.732033 4.157596 0.657 0.5115
## StatusPatient 6.061460 0.831821 7.287 1.49e-12 ***
## test_dat$Sex -0.968674 0.847673 -1.143 0.2538
## test_dat$Age 0.185190 0.213239 0.868 0.3856
## test_dat$Age2 -0.006135 0.002601 -2.359 0.0188 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.589 on 437 degrees of freedom
## Multiple R-squared: 0.3015, Adjusted R-squared: 0.2951
## F-statistic: 47.16 on 4 and 437 DF, p-value: < 2.2e-16summary(lm(Brain_PAD~status_sch+test_dat$Sex+test_dat$Age+test_dat$Age2,data=dat))##
## Call:
## lm(formula = Brain_PAD ~ status_sch + test_dat$Sex + test_dat$Age +
## test_dat$Age2, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -42.265 -5.243 0.229 5.836 23.194
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.476826 4.392507 0.792 0.4291
## status_schSchizophrenia 6.547650 0.892309 7.338 1.22e-12 ***
## test_dat$Sex -0.806285 0.900654 -0.895 0.3712
## test_dat$Age 0.127325 0.226187 0.563 0.5738
## test_dat$Age2 -0.005419 0.002730 -1.985 0.0479 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.649 on 400 degrees of freedom
## (결측으로 인하여 37개의 관측치가 삭제되었습니다.)
## Multiple R-squared: 0.2978, Adjusted R-squared: 0.2908
## F-statistic: 42.41 on 4 and 400 DF, p-value: < 2.2e-16dat_for_regression <- dat
replicated <- setdiff(which(!is.na(dat$status_trs) & !is.na(dat$status_ssd)),1:214)
dat[replicated,]## SubjID Brain_PAD Status status_ssd status_trs status_sch
## 415 T19 24.279541 Patient FE-SSD TRS Schizophrenia
## 419 T24 3.990798 Patient FE-SSD TRS Schizophreniadat_for_regression <- rbind(dat_for_regression,dat[replicated,])
dat_for_regression[replicated,]$status_ssd <- NA
dat_for_regression[(nrow(dat_for_regression)-1):nrow(dat_for_regression),]$status_trs <- NA
table(dat_for_regression$status_trs)##
## Control TRS
## 209 46table(dat_for_regression$status_ssd)##
## Control FE-SSD
## 209 72table(dat_for_regression$status_sch)##
## Control Schizophrenia
## 209 198test_dat_for_regression <- test_dat
test_dat_for_regression <- rbind(test_dat_for_regression,test_dat[replicated,])
dat_for_regression$status_ssd[is.na(dat_for_regression$status_ssd)] = dat_for_regression$status_trs[is.na(dat_for_regression$status_ssd)]
summary(lm(Brain_PAD~dat_for_regression$status_ssd+test_dat_for_regression$Sex+
test_dat_for_regression$Age+
test_dat_for_regression$Age2,data=dat_for_regression))##
## Call:
## lm(formula = Brain_PAD ~ dat_for_regression$status_ssd + test_dat_for_regression$Sex +
## test_dat_for_regression$Age + test_dat_for_regression$Age2,
## data = dat_for_regression)
##
## Residuals:
## Min 1Q Median 3Q Max
## -23.0328 -5.2602 -0.3887 5.6345 24.1575
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.579598 4.626825 1.638 0.102363
## dat_for_regression$status_ssdFE-SSD 4.295696 1.101383 3.900 0.000117 ***
## dat_for_regression$status_ssdTRS 10.094969 1.366464 7.388 1.3e-12 ***
## test_dat_for_regression$Sex -0.478314 0.912560 -0.524 0.600539
## test_dat_for_regression$Age -0.077507 0.249808 -0.310 0.756560
## test_dat_for_regression$Age2 -0.003534 0.003025 -1.169 0.243448
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.753 on 321 degrees of freedom
## (결측으로 인하여 117개의 관측치가 삭제되었습니다.)
## Multiple R-squared: 0.4011, Adjusted R-squared: 0.3917
## F-statistic: 42.99 on 5 and 321 DF, p-value: < 2.2e-16# ANOVA results (Control vs FE-SSD vs TRS)
summary(res_aov <- aov(dat_for_regression$Brain_PAD ~ dat_for_regression$status_ssd))## Df Sum Sq Mean Sq F value Pr(>F)
## dat_for_regression$status_ssd 2 4642 2321.2 27.28 1.13e-11 ***
## Residuals 324 27572 85.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 결측으로 인하여 117개의 관측치가 삭제되었습니다.pairwise.t.test(dat_for_regression$Brain_PAD,dat_for_regression$status_ssd, p.adj = "fdr")$p.value## Control FE-SSD
## FE-SSD 6.466828e-08 NA
## TRS 2.793844e-08 0.3059112pairwise.t.test(dat_for_regression$Brain_PAD,dat_for_regression$status_ssd, p.adj = "bonferroni")$p.value## Control FE-SSD
## FE-SSD 1.293366e-07 NA
## TRS 2.793844e-08 0.9177336# t.test results (Control vs Schizo)
t.test(dat$Brain_PAD~dat$status_sch)##
## Welch Two Sample t-test
##
## data: dat$Brain_PAD by dat$status_sch
## t = -6.6336, df = 385.01, p-value = 1.112e-10
## alternative hypothesis: true difference in means between group Control and group Schizophrenia is not equal to 0
## 95 percent confidence interval:
## -8.385513 -4.551188
## sample estimates:
## mean in group Control mean in group Schizophrenia
## -1.887131 4.581220t.test(dat$Brain_PAD~dat$status_sch, var.equal=T)##
## Two Sample t-test
##
## data: dat$Brain_PAD by dat$status_sch
## t = -6.6663, df = 403, p-value = 8.657e-11
## alternative hypothesis: true difference in means between group Control and group Schizophrenia is not equal to 0
## 95 percent confidence interval:
## -8.375855 -4.560847
## sample estimates:
## mean in group Control mean in group Schizophrenia
## -1.887131 4.581220#### corrected_Brain_PAD ####
dat <- data.frame(SubjID=test_dat$SubjID,corrected_Brain_PAD = c(cPAD_rid),Status=rep(c("Control","Patient"),c(nrow(Test_HC_data),nrow(Test_Patient_dat))))
dat$status_sch <- dat$status_trs <- dat$status_ssd <- rep(NA,nrow(dat))
dat$status_sch[test_dat$SubjID %in% Test_HC_data$SubjID ] <- "Control"
dat$status_sch[test_dat$SubjID %in% Test_schizophrenia_dat$SubjID ] <- "Schizophrenia"
dat$status_trs[test_dat$SubjID %in% Test_TRS_dat$SubjID ] <- "TRS"
dat$status_trs[test_dat$SubjID %in% Test_HC_data$SubjID ] <- "Control"
dat$status_ssd[test_dat$SubjID %in% Test_FESSD_dat$SubjID ] <- "FE-SSD"
dat$status_ssd[test_dat$SubjID %in% Test_HC_data$SubjID ] <- "Control"
dat %>% dplyr::select(corrected_Brain_PAD,Status) %>% group_by(Status) %>% summarise_each(funs(mean,sd,std.error))## # A tibble: 2 × 4
## Status mean sd std.error
## <chr> <dbl> <dbl> <dbl>
## 1 Control -3.09 7.57 0.524
## 2 Patient 2.77 9.39 0.615dat %>% dplyr::select(corrected_Brain_PAD,status_sch) %>% group_by(status_sch) %>% summarise_each(funs(mean,sd,std.error))## # A tibble: 3 × 4
## status_sch mean sd std.error
## <chr> <dbl> <dbl> <dbl>
## 1 Control -3.09 7.57 0.524
## 2 Schizophrenia 3.11 9.65 0.690
## 3 <NA> 0.993 7.69 1.26dat %>% dplyr::select(corrected_Brain_PAD,status_trs) %>% group_by(status_trs) %>% summarise_each(funs(mean,sd,std.error))## # A tibble: 3 × 4
## status_trs mean sd std.error
## <chr> <dbl> <dbl> <dbl>
## 1 Control -3.09 7.57 0.524
## 2 TRS 5.77 7.32 1.08
## 3 <NA> 2.04 9.70 0.710dat %>% dplyr::select(corrected_Brain_PAD,status_ssd) %>% group_by(status_ssd) %>% summarise_each(funs(mean,sd,std.error))## # A tibble: 3 × 4
## status_ssd mean sd std.error
## <chr> <dbl> <dbl> <dbl>
## 1 Control -3.09 7.57 0.524
## 2 FE-SSD 1.85 8.66 1.02
## 3 <NA> 3.19 9.69 0.764library(ggplot2)
ggplot(dat,aes(x = corrected_Brain_PAD, fill = Status)) + geom_density(alpha = 0.5) + scale_fill_manual(values=palette[c(3,5)]) +
coord_cartesian(xlim =c(-30,30)) + theme_bw()
A <- ggplot(dat[!is.na(dat$status_sch),],aes(x = corrected_Brain_PAD, fill = status_sch)) + scale_fill_manual(values=palette[c(3,4)]) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))
B <- ggplot(dat[!is.na(dat$status_ssd),],aes(x = corrected_Brain_PAD, fill = status_ssd)) + scale_fill_manual(values=palette[c(3,1)]) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))
C <- ggplot(dat[!is.na(dat$status_trs),],aes(x = corrected_Brain_PAD, fill = status_trs)) +scale_fill_manual(values=palette[c(3,2)]) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))
png("S5-2(ridge).png", width = 1400, height = 300)
ggarrange(A,B,C,
labels = c("A","B","C"),
ncol = 3, nrow = 1,
font.label = list(size=20))
dev.off()## png
## 2# ANOVA results (Control vs FE-SSD vs TRS)
dat_for_regression <- dat
replicated <- setdiff(which(!is.na(dat$status_trs) & !is.na(dat$status_ssd)),1:214)
dat[replicated,]## SubjID corrected_Brain_PAD Status status_ssd status_trs status_sch
## 981 T19 19.94609021 Patient FE-SSD TRS Schizophrenia
## 985 T24 -0.03676607 Patient FE-SSD TRS Schizophreniadat_for_regression <- rbind(dat_for_regression,dat[replicated,])
dat_for_regression[replicated,]$status_ssd <- NA
dat_for_regression[(nrow(dat_for_regression)-1):nrow(dat_for_regression),]$status_trs <- NA
table(dat_for_regression$status_trs)##
## Control TRS
## 209 46table(dat_for_regression$status_ssd)##
## Control FE-SSD
## 209 72test_dat_for_regression <- test_dat
test_dat_for_regression <- rbind(test_dat_for_regression,test_dat[replicated,])
dat_for_regression$status_ssd[is.na(dat_for_regression$status_ssd)] = dat_for_regression$status_trs[is.na(dat_for_regression$status_ssd)]
summary(res_aov <- aov(dat_for_regression$corrected_Brain_PAD~
dat_for_regression$status_ssd))## Df Sum Sq Mean Sq F value Pr(>F)
## dat_for_regression$status_ssd 2 3590 1794.8 29.58 1.59e-12 ***
## Residuals 324 19661 60.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 결측으로 인하여 117개의 관측치가 삭제되었습니다.pairwise.t.test(dat_for_regression$corrected_Brain_PAD,
dat_for_regression$status_ssd, p.adj = "fdr")$p.value## Control FE-SSD
## FE-SSD 7.577113e-06 NA
## TRS 4.792668e-11 0.007980272pairwise.t.test(dat_for_regression$corrected_Brain_PAD,
dat_for_regression$status_ssd, p.adj = "bonferroni")$p.value## Control FE-SSD
## FE-SSD 1.515423e-05 NA
## TRS 4.792668e-11 0.02394082# t.test results (Control vs Schizo)
t.test(dat$corrected_Brain_PAD~dat$status_sch)##
## Welch Two Sample t-test
##
## data: dat$corrected_Brain_PAD by dat$status_sch
## t = -7.1625, df = 369.56, p-value = 4.317e-12
## alternative hypothesis: true difference in means between group Control and group Schizophrenia is not equal to 0
## 95 percent confidence interval:
## -7.904795 -4.499344
## sample estimates:
## mean in group Control mean in group Schizophrenia
## -3.092237 3.109832t.test(dat$corrected_Brain_PAD~dat$status_sch, var.equal=T)##
## Two Sample t-test
##
## data: dat$corrected_Brain_PAD by dat$status_sch
## t = -7.2177, df = 403, p-value = 2.65e-12
## alternative hypothesis: true difference in means between group Control and group Schizophrenia is not equal to 0
## 95 percent confidence interval:
## -7.891324 -4.512815
## sample estimates:
## mean in group Control mean in group Schizophrenia
## -3.092237 3.109832summary(res_aov <- aov(dat$corrected_Brain_PAD~dat$status_sch))## Df Sum Sq Mean Sq F value Pr(>F)
## dat$status_sch 1 3891 3891 52.09 2.65e-12 ***
## Residuals 403 30098 75
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 결측으로 인하여 37개의 관측치가 삭제되었습니다.SVR[RBF] (Nested CV)
library(kernlab)##
## 다음의 패키지를 부착합니다: 'kernlab'## The following object is masked from 'package:ggplot2':
##
## alphastart.time <- Sys.time()
Nfolds <- 10
svm.fun <- function(N,i){
Nfolds <- N
set.seed(i)
outer_folds <- createFolds(Train_HC_data$Age,k=Nfolds)
cv.svm <- lapply(outer_folds,function(index) {
tr <- Train_HC_data[-index,]
ts <- Train_HC_data[index,]
svm.fit <- ksvm(as.matrix(tr[, c(3,5:81)]),
tr[,2],
kpar="automatic",
kernel="rbfdot",cross=Nfolds)
c(svm.fit@cross,svm.fit@kernelf@kpar[["sigma"]])
})
}
Ntimes <- 10
RES <- NULL
for(i in 1:Ntimes){
res <- mclapply(10,i,FUN=svm.fun,mc.cores = 6)
RES <- rbind(RES,t(data.frame(res)))
}
end.time <- Sys.time()
time.taken <- end.time - start.time
time.taken## Time difference of 29.7719 secssvm.fit.77 <- ksvm(as.matrix(Train_HC_data[,c(3,5:81)]),Train_HC_data[,2],
kpar=list(sigma=RES[,2][which.min(RES[,1])]),
kernel="rbfdot")
svm.fit.77## Support Vector Machine object of class "ksvm"
##
## SV type: eps-svr (regression)
## parameter : epsilon = 0.1 cost C = 1
##
## Gaussian Radial Basis kernel function.
## Hyperparameter : sigma = 0.0077714944308448
##
## Number of Support Vectors : 449
##
## Objective Function Value : -144.8082
## Training error : 0.141236# Train Control
ctrl_pred_age <- predict(svm.fit.77 ,as.matrix(Train_HC_data[,c(3,5:81)]))
train_error_ctrl_svm <- ctrl_pred_age-Train_HC_data$Age
par(mfrow=c(1,1))
plot(ctrl_pred_age,Train_HC_data$Age,
xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (Train control data)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",
legend=c(paste("r=",round(cor(Train_HC_data$Age,ctrl_pred_age),4)),
paste("R2=",round(cor(Train_HC_data$Age,ctrl_pred_age)^2,4)),
paste("MAE=",round(sum(abs(train_error_ctrl_svm)/length(train_error_ctrl_svm)),4)),
paste("VEcv=",round(1-sum(train_error_ctrl_svm^2/(var(Train_HC_data$Age)*(length(Train_HC_data$Age)-1))),4)*100,"%")) ,bty="n")
# Test Control
ctrl_pred_age <- predict(svm.fit.77 , as.matrix(Test_HC_data[,c(3,5:81)]))
test_error_ctrl_svm <- ctrl_pred_age-Test_HC_data$Age
plot(ctrl_pred_age,Test_HC_data$Age,
xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (Test control data)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",
legend=c(paste("r=",round(cor(Test_HC_data$Age,ctrl_pred_age),4)),
paste("R2=",round(cor(Test_HC_data$Age,ctrl_pred_age)^2,4)),
paste("MAE=",round(sum(abs(test_error_ctrl_svm)/length(test_error_ctrl_svm)),4)),
paste("VEcv=",round(1-sum(test_error_ctrl_svm^2/(var(Test_HC_data$Age)*(length(Test_HC_data$Age)-1))),4)*100,"%")) ,bty="n")
# Patient
pat_pred_age <- predict(svm.fit.77 ,as.matrix(Test_Patient_dat[,c(3,5:81)]))
test_error_pat_svm <- pat_pred_age-Test_Patient_dat$Age
plot(pat_pred_age,Test_Patient_dat$Age,xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (Patient)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",legend=c(paste("r=",round(cor(Test_Patient_dat$Age,pat_pred_age),4)),paste("R2=",round(cor(Test_Patient_dat$Age,pat_pred_age)^2,4)), paste("MAE=",round(sum(abs(test_error_pat_svm)/length(test_error_pat_svm)),4)),
paste("VEcv=",round(1-sum(test_error_pat_svm^2/(var(Test_Patient_dat$Age)*(length(Test_Patient_dat$Age)-1))),4)*100,"%")),bty="n")
# Patient - schizo
sch_pred_age <- predict(svm.fit.77, as.matrix(Test_schizophrenia_dat[,c(3,5:81)]))
test_error_sch_svm <- sch_pred_age-Test_schizophrenia_dat$Age
plot(sch_pred_age,Test_schizophrenia_dat$Age,xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (Schizophrenia)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",legend=c(paste("r=",round(cor(Test_schizophrenia_dat$Age,sch_pred_age),4)),
paste("R2=",round(cor(Test_schizophrenia_dat$Age,sch_pred_age)^2,4)),
paste("MAE=",round(sum(abs(test_error_sch_svm)/length(test_error_sch_svm)),4)),
paste("VEcv=",round(1-sum(test_error_sch_svm^2/(var(Test_schizophrenia_dat$Age)*(length(Test_schizophrenia_dat$Age)-1))),4)*100,"%")),bty="n")
# Patient - TRS
TRS_pred_age <- predict(svm.fit.77, as.matrix(Test_TRS_dat[,c(3,5:81)]))
test_error_trs_svm <- TRS_pred_age-Test_TRS_dat$Age
plot(TRS_pred_age,Test_TRS_dat$Age,xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (TRS)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",legend=c(paste("r=",round(cor(Test_TRS_dat$Age,TRS_pred_age),4)),
paste("R2=",round(cor(Test_TRS_dat$Age,TRS_pred_age)^2,4)),
paste("MAE=",round(sum(abs(test_error_trs_svm)/length(test_error_trs_svm)),4)),
paste("VEcv=",round(1-sum(test_error_trs_svm^2/(var(Test_TRS_dat$Age)*(length(Test_TRS_dat$Age)-1))),4)*100,"%")),bty="n")
# Patient - FESSD
SSD_pred_age <- predict(svm.fit.77, as.matrix(Test_FESSD_dat[,c(3,5:81)]))
test_error_ssd_svm <- SSD_pred_age-Test_FESSD_dat$Age
plot(SSD_pred_age,Test_FESSD_dat$Age,xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (FE-SSD)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",legend=c(paste("r=",round(cor(Test_FESSD_dat$Age,SSD_pred_age),4)),
paste("R2=",round(cor(Test_FESSD_dat$Age,SSD_pred_age)^2,4)),
paste("MAE=",round(sum(abs(test_error_ssd_svm)/length(test_error_ssd_svm)),4)),
paste("VEcv=",round(1-sum(test_error_ssd_svm^2/(var(Test_FESSD_dat$Age)*(length(Test_FESSD_dat$Age)-1))),4)*100,"%")),bty="n")
Brain PAD correction-plot
#### Brain PAD correction-HC
pred_age <- predict(svm.fit.77, as.matrix(Test_HC_data[,c(3,5:81)]))
test_error_svm <- pred_age-Test_HC_data$Age
PAD_svm <- test_error_svm
png("2_SVR_corr.png", width = 800, height = 600)
par(mfrow=c(2,2), oma=c(1,3,1,1))
plot(Test_HC_data$Age,PAD_svm,xlim=c(10,80),ylim=c(-40,40),pch=16,
xlab="Chronological age (Test HC)",ylab="Brain PAD uncorrected")
abline(h=0,lty=2,col="red")
abline(lm(PAD_svm~Test_HC_data$Age))
legend("topright",legend=paste("r=",round(cor(PAD_svm,Test_HC_data$Age),4),"(p<0.05)",sep=""))
cor.test(PAD_svm,Test_HC_data$Age)$p.value## [1] 3.178454e-25mtext("A", side=3, line=-3, cex=2, adj = -0.0, outer=TRUE)
Model_linear <- lm(PAD_svm~Test_HC_data$Age+Test_HC_data$Sex)
NullModel <- lm(PAD_svm~1)
anova(NullModel,Model_linear,test="Chisq") # Reject H0 (Full model is better)## Analysis of Variance Table
##
## Model 1: PAD_svm ~ 1
## Model 2: PAD_svm ~ Test_HC_data$Age + Test_HC_data$Sex
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 208 14952
## 2 206 8783 2 6169.3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Test_HC_data$Age2 <- Test_HC_data$Age^2
Model_quad <- lm(PAD_svm~+Test_HC_data$Sex+Test_HC_data$Age+Test_HC_data$Age2)
anova(NullModel,Model_quad,test="Chisq") # Reject H0 (Full model is better)## Analysis of Variance Table
##
## Model 1: PAD_svm ~ 1
## Model 2: PAD_svm ~ +Test_HC_data$Sex + Test_HC_data$Age + Test_HC_data$Age2
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 208 14952.3
## 2 205 8743.2 3 6209.1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(Model_linear,Model_quad,test="Chisq") # Not Reject H0 (Reduced model is better)## Analysis of Variance Table
##
## Model 1: PAD_svm ~ Test_HC_data$Age + Test_HC_data$Sex
## Model 2: PAD_svm ~ +Test_HC_data$Sex + Test_HC_data$Age + Test_HC_data$Age2
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 206 8783.0
## 2 205 8743.2 1 39.742 0.3344cPAD_svm <- lm(PAD_svm~+Test_HC_data$Sex+Test_HC_data$Age)$residuals
plot(Test_HC_data$Age,cPAD_svm,xlim=c(10,75),ylim=c(-40,40),pch=16,
xlab="Chronological age (Test HC)",ylab="Brain PAD corrected for age")
abline(lm(cPAD_svm~Test_HC_data$Age))
abline(h=0,lty=2,col="red")
mtext("B", side=3, line=-3, cex=2, adj = 0.5, outer=TRUE)
#### Brain PAD correction-Patient
pred_age <- predict(svm.fit.77, as.matrix(Test_Patient_dat[,c(3,5:81)]))
test_error_svm <- pred_age-Test_Patient_dat$Age
PAD_svm <- test_error_svm
plot(Test_Patient_dat$Age,PAD_svm,xlim=c(10,80),ylim=c(-40,40),pch=16,
xlab="Chronological age (Test Patient)",ylab="Brain PAD uncorrected")
abline(h=0,lty=2,col="red")
abline(lm(PAD_svm~Test_Patient_dat$Age))
legend("topright",legend=paste("r=",round(cor(PAD_svm,Test_Patient_dat$Age),4),"(p<0.05)",sep=""))
cor.test(PAD_svm,Test_Patient_dat$Age)$p.value## [1] 4.532322e-23mtext("C", side=3, line=-27, cex=2, adj = -0.0, outer=TRUE)
Model_linear <- lm(PAD_svm~Test_Patient_dat$Age+Test_Patient_dat$Sex)
NullModel <- lm(PAD_svm~1)
anova(NullModel,Model_linear,test="Chisq") # Reject H0 (Full model is better)## Analysis of Variance Table
##
## Model 1: PAD_svm ~ 1
## Model 2: PAD_svm ~ Test_Patient_dat$Age + Test_Patient_dat$Sex
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 232 23999
## 2 230 15552 2 8447.1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Test_Patient_dat$Age2 <- Test_Patient_dat$Age^2
Model_quad <- lm(PAD_svm~+Test_Patient_dat$Sex+Test_Patient_dat$Age+Test_Patient_dat$Age2)
anova(NullModel,Model_quad,test="Chisq") # Reject H0 (Full model is better)## Analysis of Variance Table
##
## Model 1: PAD_svm ~ 1
## Model 2: PAD_svm ~ +Test_Patient_dat$Sex + Test_Patient_dat$Age + Test_Patient_dat$Age2
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 232 23999
## 2 229 14610 3 9389 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(Model_linear,Model_quad,test="Chisq") # Reject H0 (Full model is better)## Analysis of Variance Table
##
## Model 1: PAD_svm ~ Test_Patient_dat$Age + Test_Patient_dat$Sex
## Model 2: PAD_svm ~ +Test_Patient_dat$Sex + Test_Patient_dat$Age + Test_Patient_dat$Age2
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 230 15552
## 2 229 14610 1 941.93 0.0001218 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cPAD_svm <- lm(PAD_svm~+Test_Patient_dat$Sex+Test_Patient_dat$Age+Test_Patient_dat$Age2)$residuals
plot(Test_Patient_dat$Age,cPAD_svm,xlim=c(10,75),ylim=c(-40,40),pch=16,
xlab="Chronological age (Test Patienit)",ylab="Brain PAD corrected for age")
abline(lm(cPAD_svm~Test_Patient_dat$Age+Test_Patient_dat$Age2))## Warning in abline(lm(cPAD_svm ~ Test_Patient_dat$Age + Test_Patient_dat$Age2)):
## only using the first two of 3 regression coefficientsabline(h=0,lty=2,col="red")
mtext("D", side=3, line=-27, cex=2, adj = 0.5, outer=TRUE)
dev.off()## png
## 2Brain PAD correction -table
test_dat <- rbind(Test_HC_data,Test_Patient_dat)
test_dat$Dx <- c(rep(0,nrow(Test_HC_data)),rep(1,nrow(Test_Patient_dat)))
test_dat$Dx <- factor(test_dat$Dx)
pred_age <- predict(svm.fit.77 , as.matrix(test_dat[,c(3,5:81)]))
test_error_svm <- pred_age-test_dat$Age
PAD_svm <- test_error_svm
mean(PAD_svm) ; sd(PAD_svm) ; median(PAD_svm)## [1] 0.07811399## [1] 9.6726## [1] 0.6968941Model_linear <- lm(PAD_svm~test_dat$Age+test_dat$Sex)
NullModel <- lm(PAD_svm~1)
anova(NullModel,Model_linear,test="Chisq") # Reject H0 (Full model is better)## Analysis of Variance Table
##
## Model 1: PAD_svm ~ 1
## Model 2: PAD_svm ~ test_dat$Age + test_dat$Sex
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 441 41260
## 2 439 26917 2 14342 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1test_dat$Age2 <- test_dat$Age^2
Model_quad <- lm(PAD_svm~+test_dat$Sex+test_dat$Age+test_dat$Age2)
anova(NullModel,Model_quad,test="Chisq") # Reject H0 (Full model is better)## Analysis of Variance Table
##
## Model 1: PAD_svm ~ 1
## Model 2: PAD_svm ~ +test_dat$Sex + test_dat$Age + test_dat$Age2
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 441 41260
## 2 438 26030 3 15230 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(Model_linear,Model_quad,test="Chisq") # Not Reject H0 (Reduced model is better)## Analysis of Variance Table
##
## Model 1: PAD_svm ~ test_dat$Age + test_dat$Sex
## Model 2: PAD_svm ~ +test_dat$Sex + test_dat$Age + test_dat$Age2
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 439 26917
## 2 438 26030 1 887.43 0.0001114 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cPAD_svm <- lm(PAD_svm~+test_dat$Sex+test_dat$Age+test_dat$Age2)$residuals
mean(cPAD_svm) ; sd(cPAD_svm) ; median(cPAD_svm)## [1] 2.730541e-16## [1] 7.682741## [1] -0.2859405library(plotrix)
dat <- data.frame(SubjID=test_dat$SubjID,Brain_PAD = c(PAD_svm),Status=rep(c("Control","Patient"),c(nrow(Test_HC_data),nrow(Test_Patient_dat))))
dat$status_sch <- dat$status_trs <- dat$status_ssd <- rep(NA,nrow(dat))
dat$status_sch[test_dat$SubjID %in% Test_HC_data$SubjID ] <- "Control"
dat$status_sch[test_dat$SubjID %in% Test_schizophrenia_dat$SubjID ] <- "Schizophrenia"
dat$status_trs[test_dat$SubjID %in% Test_TRS_dat$SubjID ] <- "TRS"
dat$status_trs[test_dat$SubjID %in% Test_HC_data$SubjID ] <- "Control"
dat$status_ssd[test_dat$SubjID %in% Test_FESSD_dat$SubjID ] <- "FE-SSD"
dat$status_ssd[test_dat$SubjID %in% Test_HC_data$SubjID ] <- "Control"
library(dplyr)
palette <- c("FE-SSD"="darkgreen","TRS"="red","Control"="grey","Schizophrenia"="blue","Patient"="orange")
dat %>% dplyr::select(Brain_PAD,Status) %>% group_by(Status) %>% summarise_each(funs(mean,sd,std.error))## # A tibble: 2 × 4
## Status mean sd std.error
## <chr> <dbl> <dbl> <dbl>
## 1 Control -2.33 8.48 0.586
## 2 Patient 2.24 10.2 0.666dat %>% dplyr::select(Brain_PAD,status_sch) %>% group_by(status_sch) %>% summarise_each(funs(mean,sd,std.error))## # A tibble: 3 × 4
## status_sch mean sd std.error
## <chr> <dbl> <dbl> <dbl>
## 1 Control -2.33 8.48 0.586
## 2 Schizophrenia 1.76 10.2 0.728
## 3 <NA> 4.79 9.83 1.62dat %>% dplyr::select(Brain_PAD,status_trs) %>% group_by(status_trs) %>% summarise_each(funs(mean,sd,std.error))## # A tibble: 3 × 4
## status_trs mean sd std.error
## <chr> <dbl> <dbl> <dbl>
## 1 Control -2.33 8.48 0.586
## 2 TRS 3.06 8.11 1.20
## 3 <NA> 2.04 10.6 0.777dat %>% dplyr::select(Brain_PAD,status_ssd) %>% group_by(status_ssd) %>% summarise_each(funs(mean,sd,std.error))## # A tibble: 3 × 4
## status_ssd mean sd std.error
## <chr> <dbl> <dbl> <dbl>
## 1 Control -2.33 8.48 0.586
## 2 FE-SSD 4.53 10.8 1.27
## 3 <NA> 1.22 9.75 0.769ggplot(dat,aes(x = Brain_PAD, fill = Status)) + geom_density(alpha = 0.5) + scale_fill_manual(values=palette) +
coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))
ggplot(dat[!is.na(dat$status_sch),],aes(x = Brain_PAD, fill = status_sch)) + scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))
ggplot(dat[!is.na(dat$status_ssd),],aes(x = Brain_PAD, fill = status_ssd)) + scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))
ggplot(dat[!is.na(dat$status_trs),],aes(x = Brain_PAD, fill = status_trs)) +scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))
summary(lm(Brain_PAD~Status+test_dat$Sex+test_dat$Age+test_dat$Age2,data=dat))##
## Call:
## lm(formula = Brain_PAD ~ Status + test_dat$Sex + test_dat$Age +
## test_dat$Age2, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.9510 -4.2134 -0.1874 5.0155 26.1833
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0002983 3.5769235 0.000 0.999933
## StatusPatient 4.5084597 0.7156445 6.300 7.28e-10 ***
## test_dat$Sex 1.6882406 0.7292822 2.315 0.021079 *
## test_dat$Age 0.1881796 0.1834570 1.026 0.305581
## test_dat$Age2 -0.0074606 0.0022376 -3.334 0.000928 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.39 on 437 degrees of freedom
## Multiple R-squared: 0.4216, Adjusted R-squared: 0.4164
## F-statistic: 79.65 on 4 and 437 DF, p-value: < 2.2e-16summary(lm(Brain_PAD~status_sch+test_dat$Sex+test_dat$Age+test_dat$Age2,data=dat))##
## Call:
## lm(formula = Brain_PAD ~ status_sch + test_dat$Sex + test_dat$Age +
## test_dat$Age2, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.4697 -4.1100 -0.2096 4.9655 25.9779
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.556947 3.737911 0.149 0.88163
## status_schSchizophrenia 4.728114 0.759332 6.227 1.21e-09 ***
## test_dat$Sex 1.738308 0.766434 2.268 0.02386 *
## test_dat$Age 0.148296 0.192479 0.770 0.44149
## test_dat$Age2 -0.006932 0.002323 -2.984 0.00302 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.36 on 400 degrees of freedom
## (결측으로 인하여 37개의 관측치가 삭제되었습니다.)
## Multiple R-squared: 0.4126, Adjusted R-squared: 0.4067
## F-statistic: 70.24 on 4 and 400 DF, p-value: < 2.2e-16dat_for_regression <- dat
replicated <- setdiff(which(!is.na(dat$status_trs) & !is.na(dat$status_ssd)),1:214)
dat[replicated,]## SubjID Brain_PAD Status status_ssd status_trs status_sch
## 415 T19 22.432483 Patient FE-SSD TRS Schizophrenia
## 419 T24 4.575506 Patient FE-SSD TRS Schizophreniadat_for_regression <- rbind(dat_for_regression,dat[replicated,])
dat_for_regression[replicated,]$status_ssd <- NA
dat_for_regression[(nrow(dat_for_regression)-1):nrow(dat_for_regression),]$status_trs <- NA
table(dat_for_regression$status_trs)##
## Control TRS
## 209 46table(dat_for_regression$status_ssd)##
## Control FE-SSD
## 209 72test_dat_for_regression <- test_dat
test_dat_for_regression <- rbind(test_dat_for_regression,test_dat[replicated,])
dat_for_regression$status_ssd[is.na(dat_for_regression$status_ssd)] = dat_for_regression$status_trs[is.na(dat_for_regression$status_ssd)]
summary(lm(Brain_PAD~dat_for_regression$status_ssd+test_dat_for_regression$Sex+
test_dat_for_regression$Age+
test_dat_for_regression$Age2,data=dat_for_regression))##
## Call:
## lm(formula = Brain_PAD ~ dat_for_regression$status_ssd + test_dat_for_regression$Sex +
## test_dat_for_regression$Age + test_dat_for_regression$Age2,
## data = dat_for_regression)
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.1642 -4.2029 -0.3043 4.9129 20.1372
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.052667 4.128789 1.224 0.221939
## dat_for_regression$status_ssdFE-SSD 3.628643 0.982829 3.692 0.000261 ***
## dat_for_regression$status_ssdTRS 7.591653 1.219377 6.226 1.5e-09 ***
## test_dat_for_regression$Sex 1.954197 0.814331 2.400 0.016975 *
## test_dat_for_regression$Age -0.110802 0.222918 -0.497 0.619492
## test_dat_for_regression$Age2 -0.003924 0.002699 -1.454 0.146950
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.918 on 321 degrees of freedom
## (결측으로 인하여 117개의 관측치가 삭제되었습니다.)
## Multiple R-squared: 0.4736, Adjusted R-squared: 0.4654
## F-statistic: 57.77 on 5 and 321 DF, p-value: < 2.2e-16# ANOVA results (Control vs FE-SSD vs TRS)
summary(res_aov <- aov(dat_for_regression$Brain_PAD ~ dat_for_regression$status_ssd))## Df Sum Sq Mean Sq F value Pr(>F)
## dat_for_regression$status_ssd 2 3043 1521.5 18.85 1.8e-08 ***
## Residuals 324 26146 80.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 결측으로 인하여 117개의 관측치가 삭제되었습니다.pairwise.t.test(dat_for_regression$Brain_PAD,dat_for_regression$status_ssd, p.adj = "fdr")$p.value## Control FE-SSD
## FE-SSD 1.453155e-07 NA
## TRS 4.012511e-04 0.3868526pairwise.t.test(dat_for_regression$Brain_PAD,dat_for_regression$status_ssd, p.adj = "bonferroni")$p.value## Control FE-SSD
## FE-SSD 1.453155e-07 NA
## TRS 8.025021e-04 1# t.test results (Control vs Schizo)
t.test(dat$Brain_PAD~dat$status_sch)##
## Welch Two Sample t-test
##
## data: dat$Brain_PAD by dat$status_sch
## t = -4.383, df = 380.21, p-value = 1.516e-05
## alternative hypothesis: true difference in means between group Control and group Schizophrenia is not equal to 0
## 95 percent confidence interval:
## -5.933928 -2.258724
## sample estimates:
## mean in group Control mean in group Schizophrenia
## -2.334937 1.761390t.test(dat$Brain_PAD~dat$status_sch, var.equal=T)##
## Two Sample t-test
##
## data: dat$Brain_PAD by dat$status_sch
## t = -4.4087, df = 403, p-value = 1.335e-05
## alternative hypothesis: true difference in means between group Control and group Schizophrenia is not equal to 0
## 95 percent confidence interval:
## -5.922893 -2.269759
## sample estimates:
## mean in group Control mean in group Schizophrenia
## -2.334937 1.761390#### corrected_Brain_PAD ####
dat <- data.frame(SubjID=test_dat$SubjID,corrected_Brain_PAD = c(cPAD_svm),Status=rep(c("Control","Patient"),c(nrow(Test_HC_data),nrow(Test_Patient_dat))))
dat$status_sch <- dat$status_trs <- dat$status_ssd <- rep(NA,nrow(dat))
dat$status_sch[test_dat$SubjID %in% Test_HC_data$SubjID ] <- "Control"
dat$status_sch[test_dat$SubjID %in% Test_schizophrenia_dat$SubjID ] <- "Schizophrenia"
dat$status_trs[test_dat$SubjID %in% Test_TRS_dat$SubjID ] <- "TRS"
dat$status_trs[test_dat$SubjID %in% Test_HC_data$SubjID ] <- "Control"
dat$status_ssd[test_dat$SubjID %in% Test_FESSD_dat$SubjID ] <- "FE-SSD"
dat$status_ssd[test_dat$SubjID %in% Test_HC_data$SubjID ] <- "Control"
dat %>% dplyr::select(corrected_Brain_PAD,Status) %>% group_by(Status) %>% summarise_each(funs(mean,sd,std.error))## # A tibble: 2 × 4
## Status mean sd std.error
## <chr> <dbl> <dbl> <dbl>
## 1 Control -2.30 6.59 0.456
## 2 Patient 2.06 8.02 0.525dat %>% dplyr::select(corrected_Brain_PAD,status_sch) %>% group_by(status_sch) %>% summarise_each(funs(mean,sd,std.error))## # A tibble: 3 × 4
## status_sch mean sd std.error
## <chr> <dbl> <dbl> <dbl>
## 1 Control -2.30 6.59 0.456
## 2 Schizophrenia 2.19 8.08 0.577
## 3 <NA> 1.40 7.74 1.27dat %>% dplyr::select(corrected_Brain_PAD,status_trs) %>% group_by(status_trs) %>% summarise_each(funs(mean,sd,std.error))## # A tibble: 3 × 4
## status_trs mean sd std.error
## <chr> <dbl> <dbl> <dbl>
## 1 Control -2.30 6.59 0.456
## 2 TRS 4.34 6.32 0.932
## 3 <NA> 1.50 8.30 0.607dat %>% dplyr::select(corrected_Brain_PAD,status_ssd) %>% group_by(status_ssd) %>% summarise_each(funs(mean,sd,std.error))## # A tibble: 3 × 4
## status_ssd mean sd std.error
## <chr> <dbl> <dbl> <dbl>
## 1 Control -2.30 6.59 0.456
## 2 FE-SSD 1.67 8.18 0.964
## 3 <NA> 2.24 7.96 0.627ggplot(dat,aes(x = corrected_Brain_PAD, fill = Status)) + geom_density(alpha = 0.5) + scale_fill_manual(values=palette) +
coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))
ggplot(dat[!is.na(dat$status_sch),],aes(x = corrected_Brain_PAD, fill = status_sch)) + scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))
ggplot(dat[!is.na(dat$status_ssd),],aes(x = corrected_Brain_PAD, fill = status_ssd)) + scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))
ggplot(dat[!is.na(dat$status_trs),],aes(x = corrected_Brain_PAD, fill = status_trs)) +scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))
library(ggplot2)
ggplot(dat,aes(x = corrected_Brain_PAD, fill = Status)) + geom_density(alpha = 0.5) + scale_fill_manual(values=palette[c(3,5)]) +
coord_cartesian(xlim =c(-30,30)) + theme_bw()
A <- ggplot(dat[!is.na(dat$status_sch),],aes(x = corrected_Brain_PAD, fill = status_sch)) + scale_fill_manual(values=palette[c(3,4)]) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))
B <- ggplot(dat[!is.na(dat$status_ssd),],aes(x = corrected_Brain_PAD, fill = status_ssd)) + scale_fill_manual(values=palette[c(3,1)]) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))
C <- ggplot(dat[!is.na(dat$status_trs),],aes(x = corrected_Brain_PAD, fill = status_trs)) +scale_fill_manual(values=palette[c(3,2)]) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))
png("p2(SVR).png", width = 1400, height = 300)
ggarrange(A,B,C,
labels = c("A","B","C"),
ncol = 3, nrow = 1,
font.label = list(size=20))
dev.off()## png
## 2dat_for_regression <- dat
replicated <- setdiff(which(!is.na(dat$status_trs) & !is.na(dat$status_ssd)),1:214)
dat[replicated,]## SubjID corrected_Brain_PAD Status status_ssd status_trs status_sch
## 415 T19 16.7968216 Patient FE-SSD TRS Schizophrenia
## 419 T24 -0.3460116 Patient FE-SSD TRS Schizophreniadat_for_regression <- rbind(dat_for_regression,dat[replicated,])
dat_for_regression[replicated,]$status_ssd <- NA
dat_for_regression[(nrow(dat_for_regression)-1):nrow(dat_for_regression),]$status_trs <- NA
table(dat_for_regression$status_trs)##
## Control TRS
## 209 46table(dat_for_regression$status_ssd)##
## Control FE-SSD
## 209 72test_dat_for_regression <- test_dat
test_dat_for_regression <- rbind(test_dat_for_regression,test_dat[replicated,])
dat_for_regression$status_ssd[is.na(dat_for_regression$status_ssd)] = dat_for_regression$status_trs[is.na(dat_for_regression$status_ssd)]
summary(res_aov <- aov(dat_for_regression$corrected_Brain_PAD~
dat_for_regression$status_ssd))## Df Sum Sq Mean Sq F value Pr(>F)
## dat_for_regression$status_ssd 2 2094 1046.9 21.77 1.34e-09 ***
## Residuals 324 15577 48.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 결측으로 인하여 117개의 관측치가 삭제되었습니다.pairwise.t.test(dat_for_regression$corrected_Brain_PAD,
dat_for_regression$status_ssd, p.adj = "fdr")$p.value## Control FE-SSD
## FE-SSD 5.380789e-05 NA
## TRS 3.071392e-08 0.04227491pairwise.t.test(dat_for_regression$corrected_Brain_PAD,
dat_for_regression$status_ssd, p.adj = "bonferroni")$p.value## Control FE-SSD
## FE-SSD 1.076158e-04 NA
## TRS 3.071392e-08 0.1268247# t.test results (Control vs Schizo)
t.test(dat$corrected_Brain_PAD~dat$status_sch)##
## Welch Two Sample t-test
##
## data: dat$corrected_Brain_PAD by dat$status_sch
## t = -6.1046, df = 376.68, p-value = 2.565e-09
## alternative hypothesis: true difference in means between group Control and group Schizophrenia is not equal to 0
## 95 percent confidence interval:
## -5.934865 -3.043082
## sample estimates:
## mean in group Control mean in group Schizophrenia
## -2.299979 2.188995t.test(dat$corrected_Brain_PAD~dat$status_sch, var.equal=T)##
## Two Sample t-test
##
## data: dat$corrected_Brain_PAD by dat$status_sch
## t = -6.1443, df = 403, p-value = 1.93e-09
## alternative hypothesis: true difference in means between group Control and group Schizophrenia is not equal to 0
## 95 percent confidence interval:
## -5.925225 -3.052722
## sample estimates:
## mean in group Control mean in group Schizophrenia
## -2.299979 2.188995summary(res_aov <- aov(dat$corrected_Brain_PAD~dat$status_sch))## Df Sum Sq Mean Sq F value Pr(>F)
## dat$status_sch 1 2038 2038 37.75 1.93e-09 ***
## Residuals 403 21757 54
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 결측으로 인하여 37개의 관측치가 삭제되었습니다.RVR[RBF] (Nested CV)
start.time <- Sys.time()
Nfolds <- 10
rvm.fun <- function(N,i){
Nfolds <- N
set.seed(i)
outer_folds <- createFolds(Train_HC_data$Age,k=Nfolds)
cv.rvm <- lapply(outer_folds,function(index) {
tr <- Train_HC_data[-index,]
ts <- Train_HC_data[index,]
rvm.fit <- rvm(as.matrix(tr[, c(3,5:81)]),
tr[,2],type="regression",
kpar="automatic",
kernel="rbfdot",cross=Nfolds)
c(rvm.fit@cross,rvm.fit@kernelf@kpar[["sigma"]])
})
}
Ntimes <- 10
RES <- NULL
for(i in 1:Ntimes){
res <- mclapply(10,i,FUN=rvm.fun,mc.cores = 6)
RES <- rbind(RES,t(data.frame(res)))
}## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
## Using automatic sigma estimation (sigest) for RBF or laplace kernelend.time <- Sys.time()
time.taken <- end.time - start.time
time.taken## Time difference of 29.99672 minsrvm.fit.77 <- rvm(as.matrix(Train_HC_data[,c(3,5:81)]),Train_HC_data[,2],
kpar=list(sigma=RES[,2][which.min(RES[,1])]),
type="regression",kernel="rbfdot")
rvm.fit.77## Relevance Vector Machine object of class "rvm"
## Problem type: regression
##
## Gaussian Radial Basis kernel function.
## Hyperparameter : sigma = 2.60538261599689e-10
##
## Number of Relevance Vectors : 33
## Variance : 105.0113
## Training error : 99.079268256# Train Control
ctrl_pred_age <- predict(rvm.fit.77 ,as.matrix(Train_HC_data[,c(3,5:81)]))
train_error_ctrl_rvm <- ctrl_pred_age-Train_HC_data$Age
par(mfrow=c(1,1))
plot(ctrl_pred_age,Train_HC_data$Age,
xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (Train control data)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",
legend=c(paste("r=",round(cor(Train_HC_data$Age,ctrl_pred_age),4)),
paste("R2=",round(cor(Train_HC_data$Age,ctrl_pred_age)^2,4)),
paste("MAE=",round(sum(abs(train_error_ctrl_rid)/length(train_error_ctrl_rid)),4)),
paste("VEcv=",round(1-sum(train_error_ctrl_rid^2/(var(Train_HC_data$Age)*(length(Train_HC_data$Age)-1))),4)*100,"%")) ,bty="n")
# Test Control
ctrl_pred_age <- predict(rvm.fit.77 , as.matrix(Test_HC_data[,c(3,5:81)]))
test_error_ctrl_rvm <- ctrl_pred_age-Test_HC_data$Age
plot(ctrl_pred_age,Test_HC_data$Age,
xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (Test control data)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",
legend=c(paste("r=",round(cor(Test_HC_data$Age,ctrl_pred_age),4)),
paste("R2=",round(cor(Test_HC_data$Age,ctrl_pred_age)^2,4)),
paste("MAE=",round(sum(abs(test_error_ctrl_rvm)/length(test_error_ctrl_rvm)),4)),
paste("VEcv=",round(1-sum(test_error_ctrl_rvm^2/(var(Test_HC_data$Age)*(length(Test_HC_data$Age)-1))),4)*100,"%")) ,bty="n")
# Patient
pat_pred_age <- predict(rvm.fit.77 ,as.matrix(Test_Patient_dat[,c(3,5:81)]))
test_error_pat_rvm <- pat_pred_age-Test_Patient_dat$Age
plot(pat_pred_age,Test_Patient_dat$Age,xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (Patient)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",legend=c(paste("r=",round(cor(Test_Patient_dat$Age,pat_pred_age),4)),paste("R2=",round(cor(Test_Patient_dat$Age,pat_pred_age)^2,4)),
paste("MAE=",round(sum(abs(test_error_pat_rid)/length(test_error_pat_rid)),4)),
paste("VEcv=",round(1-sum(test_error_pat_rid^2/(var(Test_Patient_dat$Age)*(length(Test_Patient_dat$Age)-1))),4)*100,"%")),bty="n")
# Patient - schizo
sch_pred_age <- predict(rvm.fit.77, as.matrix(Test_schizophrenia_dat[,c(3,5:81)]))
test_error_sch_rvm <- sch_pred_age-Test_schizophrenia_dat$Age
plot(sch_pred_age,Test_schizophrenia_dat$Age,xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (Schizophrenia)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",legend=c(paste("r=",round(cor(Test_schizophrenia_dat$Age,sch_pred_age),4)),
paste("R2=",round(cor(Test_schizophrenia_dat$Age,sch_pred_age)^2,4)),
paste("MAE=",round(sum(abs(test_error_sch_rvm)/length(test_error_sch_rvm)),4)),
paste("VEcv=",round(1-sum(test_error_sch_rvm^2/(var(Test_schizophrenia_dat$Age)*(length(Test_schizophrenia_dat$Age)-1))),4)*100,"%")),bty="n")
# Patient - TRS
TRS_pred_age <- predict(rvm.fit.77, as.matrix(Test_TRS_dat[,c(3,5:81)]))
test_error_trs_rvm <- TRS_pred_age-Test_TRS_dat$Age
plot(TRS_pred_age,Test_TRS_dat$Age,xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (TRS)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",legend=c(paste("r=",round(cor(Test_TRS_dat$Age,TRS_pred_age),4)),
paste("R2=",round(cor(Test_TRS_dat$Age,TRS_pred_age)^2,4)),
paste("MAE=",round(sum(abs(test_error_trs_rvm)/length(test_error_trs_rvm)),4)),
paste("VEcv=",round(1-sum(test_error_trs_rvm^2/(var(Test_TRS_dat$Age)*(length(Test_TRS_dat$Age)-1))),4)*100,"%")),bty="n")
# Patient - FESSD
SSD_pred_age <- predict(rvm.fit.77, as.matrix(Test_FESSD_dat[,c(3,5:81)]))
test_error_ssd_rvm <- SSD_pred_age-Test_FESSD_dat$Age
plot(SSD_pred_age,Test_FESSD_dat$Age,xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (FE-SSD)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",legend=c(paste("r=",round(cor(Test_FESSD_dat$Age,SSD_pred_age),4)),
paste("R2=",round(cor(Test_FESSD_dat$Age,SSD_pred_age)^2,4)),
paste("MAE=",round(sum(abs(test_error_ssd_rvm)/length(test_error_ssd_rvm)),4)),
paste("VEcv=",round(1-sum(test_error_ssd_rvm^2/(var(Test_FESSD_dat$Age)*(length(Test_FESSD_dat$Age)-1))),4)*100,"%")),bty="n")
Brain PAD correction-plot
#### Brain PAD correction-HC
pred_age <- predict(rvm.fit.77 , as.matrix(Test_HC_data[,c(3,5:81)]))
test_error_rvm <- pred_age-Test_HC_data$Age
PAD_rvm <- test_error_rvm
png("3_RVR_corr.png", width = 800, height = 600)
par(mfrow=c(2,2), oma=c(1,3,1,1))
plot(Test_HC_data$Age,PAD_rvm,xlim=c(10,80),ylim=c(-40,40),pch=16,
xlab="Chronological age (Test HC)",ylab="Brain PAD uncorrected")
abline(h=0,lty=2,col="red")
abline(lm(PAD_rvm~Test_HC_data$Age))
legend("topright",legend=paste("r=",round(cor(PAD_rvm,Test_HC_data$Age),4),"(p<0.05)",sep=""))
cor.test(PAD_rvm,Test_HC_data$Age)$p.value## [1] 9.523509e-56mtext("A", side=3, line=-3, cex=2, adj = -0.0, outer=TRUE)
Model_linear <- lm(PAD_rvm~Test_HC_data$Age+Test_HC_data$Sex)
NullModel <- lm(PAD_rvm~1)
anova(NullModel,Model_linear,test="Chisq") # Reject H0 (Full model is better)## Analysis of Variance Table
##
## Model 1: PAD_rvm ~ 1
## Model 2: PAD_rvm ~ Test_HC_data$Age + Test_HC_data$Sex
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 208 30264.1
## 2 206 8630.2 2 21634 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Test_HC_data$Age2 <- Test_HC_data$Age^2
Model_quad <- lm(PAD_rvm~Test_HC_data$Sex+Test_HC_data$Age+Test_HC_data$Age2)
anova(NullModel,Model_quad,test="Chisq") # Reject H0 (Full model is better)## Analysis of Variance Table
##
## Model 1: PAD_rvm ~ 1
## Model 2: PAD_rvm ~ Test_HC_data$Sex + Test_HC_data$Age + Test_HC_data$Age2
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 208 30264.1
## 2 205 8451.4 3 21813 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(Model_linear,Model_quad,test="Chisq") # Reject H0 (Full model is better)## Analysis of Variance Table
##
## Model 1: PAD_rvm ~ Test_HC_data$Age + Test_HC_data$Sex
## Model 2: PAD_rvm ~ Test_HC_data$Sex + Test_HC_data$Age + Test_HC_data$Age2
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 206 8630.2
## 2 205 8451.4 1 178.74 0.03732 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# cPAD_rvm <- lm(PAD_rvm~Test_HC_data$Sex+Test_HC_data$Age)$residuals
cPAD_rvm <- lm(PAD_rvm~Test_HC_data$Sex+Test_HC_data$Age+Test_HC_data$Age2)$residuals
plot(Test_HC_data$Age,cPAD_rvm,xlim=c(10,75),ylim=c(-40,40),pch=16,
xlab="Chronological age (Test HC)",ylab="Brain PAD corrected for age")
abline(lm(cPAD_rvm~Test_HC_data$Age))
abline(h=0,lty=2,col="red")
mtext("B", side=3, line=-3, cex=2, adj = 0.5, outer=TRUE)
#### Brain PAD correction-Patient
pred_age <- predict(rvm.fit.77 , as.matrix(Test_Patient_dat[,c(3,5:81)]))
test_error_rvm <- pred_age-Test_Patient_dat$Age
PAD_rvm <- test_error_rvm
plot(Test_Patient_dat$Age,PAD_rvm,xlim=c(10,80),ylim=c(-40,40),pch=16,
xlab="Chronological age (Test Patient)",ylab="Brain PAD uncorrected")
abline(h=0,lty=2,col="red")
abline(lm(PAD_rvm~Test_Patient_dat$Age))
legend("topright",legend=paste("r=",round(cor(PAD_rvm,Test_Patient_dat$Age),4),"(p<0.05)",sep=""))
cor.test(PAD_rvm,Test_Patient_dat$Age)$p.value## [1] 7.696674e-32mtext("C", side=3, line=-27, cex=2, adj = -0.0, outer=TRUE)
Model_linear <- lm(PAD_rvm~Test_Patient_dat$Age+Test_Patient_dat$Sex)
NullModel <- lm(PAD_rvm~1)
anova(NullModel,Model_linear,test="Chisq") # Reject H0 (Full model is better)## Analysis of Variance Table
##
## Model 1: PAD_rvm ~ 1
## Model 2: PAD_rvm ~ Test_Patient_dat$Age + Test_Patient_dat$Sex
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 232 33880
## 2 230 17664 2 16216 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Test_Patient_dat$Age2 <- Test_Patient_dat$Age^2
Model_quad <- lm(PAD_rvm~Test_Patient_dat$Sex+Test_Patient_dat$Age+Test_Patient_dat$Age2)
anova(NullModel,Model_quad,test="Chisq") # Reject H0 (Full model is better)## Analysis of Variance Table
##
## Model 1: PAD_rvm ~ 1
## Model 2: PAD_rvm ~ Test_Patient_dat$Sex + Test_Patient_dat$Age + Test_Patient_dat$Age2
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 232 33880
## 2 229 17628 3 16252 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(Model_linear,Model_quad,test="Chisq") # Not Reject H0 (Reduced model is better)## Analysis of Variance Table
##
## Model 1: PAD_rvm ~ Test_Patient_dat$Age + Test_Patient_dat$Sex
## Model 2: PAD_rvm ~ Test_Patient_dat$Sex + Test_Patient_dat$Age + Test_Patient_dat$Age2
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 230 17664
## 2 229 17628 1 36.347 0.492cPAD_rvm <- lm(PAD_rvm~Test_Patient_dat$Sex+Test_Patient_dat$Age)$residuals
plot(Test_Patient_dat$Age,cPAD_rvm,xlim=c(10,75),ylim=c(-40,40),pch=16,
xlab="Chronological age (Test Patient)",ylab="Brain PAD corrected for age")
abline(lm(cPAD_rvm~Test_Patient_dat$Age))
abline(h=0,lty=2,col="red")
mtext("D", side=3, line=-27, cex=2, adj = 0.5, outer=TRUE)
dev.off()## png
## 2Brain PAD correction -table
##### Brain PAD correction
test_dat <- rbind(Test_HC_data,Test_Patient_dat)
test_dat$Dx <- c(rep(0,nrow(Test_HC_data)),rep(1,nrow(Test_Patient_dat)))
test_dat$Dx <- factor(test_dat$Dx)
pred_age <- predict(rvm.fit.77 , as.matrix(test_dat[,c(3,5:81)]))
test_error_rvm <- pred_age-test_dat$Age
PAD_rvm <- test_error_rvm
mean(PAD_rvm) ; sd(PAD_rvm) ; median(PAD_rvm)## [1] -1.815347## [1] 12.19153## [1] -0.6313292Model_linear <- lm(PAD_rvm~test_dat$Age+test_dat$Sex)
NullModel <- lm(PAD_rvm~1)
anova(NullModel,Model_linear,test="Chisq") # Reject H0 (Full model is better)## Analysis of Variance Table
##
## Model 1: PAD_rvm ~ 1
## Model 2: PAD_rvm ~ test_dat$Age + test_dat$Sex
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 441 65547
## 2 439 28222 2 37325 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1test_dat$Age2 <- test_dat$Age^2
Model_quad <- lm(PAD_rvm~test_dat$Sex+test_dat$Age+test_dat$Age2)
anova(NullModel,Model_quad,test="Chisq") # Reject H0 (Full model is better)## Analysis of Variance Table
##
## Model 1: PAD_rvm ~ 1
## Model 2: PAD_rvm ~ test_dat$Sex + test_dat$Age + test_dat$Age2
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 441 65547
## 2 438 28153 3 37394 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(Model_linear,Model_quad,test="Chisq") # Not Reject H0 (Reduced model is better)## Analysis of Variance Table
##
## Model 1: PAD_rvm ~ test_dat$Age + test_dat$Sex
## Model 2: PAD_rvm ~ test_dat$Sex + test_dat$Age + test_dat$Age2
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 439 28222
## 2 438 28153 1 68.723 0.3011cPAD_rvm <- lm(PAD_rvm~test_dat$Sex+test_dat$Age)$residuals
mean(cPAD_rvm) ; sd(cPAD_rvm) ; median(cPAD_rvm)## [1] -4.868136e-16## [1] 7.999728## [1] -0.8784845library(plotrix)
dat <- data.frame(SubjID=test_dat$SubjID,Brain_PAD = c(PAD_rvm),Status=rep(c("Control","Patient"),c(nrow(Test_HC_data),nrow(Test_Patient_dat))))
dat$status_sch <- dat$status_trs <- dat$status_ssd <- rep(NA,nrow(dat))
dat$status_sch[test_dat$SubjID %in% Test_HC_data$SubjID ] <- "Control"
dat$status_sch[test_dat$SubjID %in% Test_schizophrenia_dat$SubjID ] <- "Schizophrenia"
dat$status_trs[test_dat$SubjID %in% Test_TRS_dat$SubjID ] <- "TRS"
dat$status_trs[test_dat$SubjID %in% Test_HC_data$SubjID ] <- "Control"
dat$status_ssd[test_dat$SubjID %in% Test_FESSD_dat$SubjID ] <- "FE-SSD"
dat$status_ssd[test_dat$SubjID %in% Test_HC_data$SubjID ] <- "Control"
library(dplyr)
palette <- c("FE-SSD"="darkgreen","TRS"="red","Control"="grey","Schizophrenia"="blue","Patient"="orange")
dat %>% dplyr::select(Brain_PAD,Status) %>% group_by(Status) %>% summarise_each(funs(mean,sd,std.error))## # A tibble: 2 × 4
## Status mean sd std.error
## <chr> <dbl> <dbl> <dbl>
## 1 Control -3.70 12.1 0.834
## 2 Patient -0.128 12.1 0.792dat %>% dplyr::select(Brain_PAD,status_sch) %>% group_by(status_sch) %>% summarise_each(funs(mean,sd,std.error))## # A tibble: 3 × 4
## status_sch mean sd std.error
## <chr> <dbl> <dbl> <dbl>
## 1 Control -3.70 12.1 0.834
## 2 Schizophrenia -1.51 11.6 0.831
## 3 <NA> 7.20 11.9 1.96dat %>% dplyr::select(Brain_PAD,status_trs) %>% group_by(status_trs) %>% summarise_each(funs(mean,sd,std.error))## # A tibble: 3 × 4
## status_trs mean sd std.error
## <chr> <dbl> <dbl> <dbl>
## 1 Control -3.70 12.1 0.834
## 2 TRS -2.65 11.4 1.69
## 3 <NA> 0.493 12.2 0.891dat %>% dplyr::select(Brain_PAD,status_ssd) %>% group_by(status_ssd) %>% summarise_each(funs(mean,sd,std.error))## # A tibble: 3 × 4
## status_ssd mean sd std.error
## <chr> <dbl> <dbl> <dbl>
## 1 Control -3.70 12.1 0.834
## 2 FE-SSD 5.39 12.3 1.45
## 3 <NA> -2.59 11.2 0.881ggplot(dat,aes(x = Brain_PAD, fill = Status)) + geom_density(alpha = 0.5) + scale_fill_manual(values=palette) +
coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))
ggplot(dat[!is.na(dat$status_sch),],aes(x = Brain_PAD, fill = status_sch)) + scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))
ggplot(dat[!is.na(dat$status_ssd),],aes(x = Brain_PAD, fill = status_ssd)) + scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))
ggplot(dat[!is.na(dat$status_trs),],aes(x = Brain_PAD, fill = status_trs)) +scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))
summary(lm(Brain_PAD~Status+test_dat$Sex+test_dat$Age,data=dat))##
## Call:
## lm(formula = Brain_PAD ~ Status + test_dat$Sex + test_dat$Age,
## data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -25.113 -5.046 -0.857 4.577 33.873
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.61240 1.55032 10.070 < 2e-16 ***
## StatusPatient 4.04615 0.74835 5.407 1.06e-07 ***
## test_dat$Sex 3.73670 0.76567 4.880 1.49e-06 ***
## test_dat$Age -0.67469 0.02701 -24.978 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.772 on 438 degrees of freedom
## Multiple R-squared: 0.5964, Adjusted R-squared: 0.5936
## F-statistic: 215.7 on 3 and 438 DF, p-value: < 2.2e-16summary(lm(Brain_PAD~status_sch+test_dat$Sex+test_dat$Age,data=dat))##
## Call:
## lm(formula = Brain_PAD ~ status_sch + test_dat$Sex + test_dat$Age,
## data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -25.243 -5.090 -0.875 4.716 33.706
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.95170 1.59412 9.379 < 2e-16 ***
## status_schSchizophrenia 4.03490 0.78702 5.127 4.59e-07 ***
## test_dat$Sex 3.96796 0.80572 4.925 1.24e-06 ***
## test_dat$Age -0.66716 0.02858 -23.342 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.74 on 401 degrees of freedom
## (결측으로 인하여 37개의 관측치가 삭제되었습니다.)
## Multiple R-squared: 0.5797, Adjusted R-squared: 0.5765
## F-statistic: 184.3 on 3 and 401 DF, p-value: < 2.2e-16dat_for_regression <- dat
replicated <- setdiff(which(!is.na(dat$status_trs) & !is.na(dat$status_ssd)),1:214)
dat[replicated,]## SubjID Brain_PAD Status status_ssd status_trs status_sch
## 415 T19 12.636224 Patient FE-SSD TRS Schizophrenia
## 419 T24 4.598935 Patient FE-SSD TRS Schizophreniadat_for_regression <- rbind(dat_for_regression,dat[replicated,])
dat_for_regression[replicated,]$status_ssd <- NA
dat_for_regression[(nrow(dat_for_regression)-1):nrow(dat_for_regression),]$status_trs <- NA
table(dat_for_regression$status_trs)##
## Control TRS
## 209 46table(dat_for_regression$status_ssd)##
## Control FE-SSD
## 209 72test_dat_for_regression <- test_dat
test_dat_for_regression <- rbind(test_dat_for_regression,test_dat[replicated,])
dat_for_regression$status_ssd[is.na(dat_for_regression$status_ssd)] = dat_for_regression$status_trs[is.na(dat_for_regression$status_ssd)]
summary(lm(Brain_PAD~dat_for_regression$status_ssd+test_dat_for_regression$Sex+
test_dat_for_regression$Age,data=dat_for_regression))##
## Call:
## lm(formula = Brain_PAD ~ dat_for_regression$status_ssd + test_dat_for_regression$Sex +
## test_dat_for_regression$Age, data = dat_for_regression)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.0574 -4.7342 -0.7179 4.3585 28.0182
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 18.55111 1.65831 11.187 < 2e-16 ***
## dat_for_regression$status_ssdFE-SSD 3.30188 1.00980 3.270 0.001192 **
## dat_for_regression$status_ssdTRS 5.53278 1.20911 4.576 6.78e-06 ***
## test_dat_for_regression$Sex 3.13507 0.84218 3.723 0.000233 ***
## test_dat_for_regression$Age -0.72666 0.03006 -24.177 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.18 on 322 degrees of freedom
## (결측으로 인하여 117개의 관측치가 삭제되었습니다.)
## Multiple R-squared: 0.6766, Adjusted R-squared: 0.6726
## F-statistic: 168.4 on 4 and 322 DF, p-value: < 2.2e-16# ANOVA results (Control vs FE-SSD vs TRS)
summary(res_aov <- aov(dat_for_regression$Brain_PAD ~ dat_for_regression$status_ssd))## Df Sum Sq Mean Sq F value Pr(>F)
## dat_for_regression$status_ssd 2 4483 2241.7 15.5 3.72e-07 ***
## Residuals 324 46850 144.6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 결측으로 인하여 117개의 관측치가 삭제되었습니다.pairwise.t.test(dat_for_regression$Brain_PAD,dat_for_regression$status_ssd, p.adj = "fdr")$p.value## Control FE-SSD
## FE-SSD 2.004933e-07 NA
## TRS 5.943715e-01 0.0006832325pairwise.t.test(dat_for_regression$Brain_PAD,dat_for_regression$status_ssd, p.adj = "bonferroni")$p.value## Control FE-SSD
## FE-SSD 2.004933e-07 NA
## TRS 1.000000e+00 0.001366465# t.test results (Control vs Schizo)
t.test(dat$Brain_PAD~dat$status_sch)##
## Welch Two Sample t-test
##
## data: dat$Brain_PAD by dat$status_sch
## t = -1.8556, df = 402.68, p-value = 0.06425
## alternative hypothesis: true difference in means between group Control and group Schizophrenia is not equal to 0
## 95 percent confidence interval:
## -4.5005093 0.1299218
## sample estimates:
## mean in group Control mean in group Schizophrenia
## -3.696547 -1.511253t.test(dat$Brain_PAD~dat$status_sch, var.equal=T)##
## Two Sample t-test
##
## data: dat$Brain_PAD by dat$status_sch
## t = -1.8534, df = 403, p-value = 0.06455
## alternative hypothesis: true difference in means between group Control and group Schizophrenia is not equal to 0
## 95 percent confidence interval:
## -4.5031844 0.1325969
## sample estimates:
## mean in group Control mean in group Schizophrenia
## -3.696547 -1.511253#### corrected_Brain_PAD ####
dat <- data.frame(SubjID=test_dat$SubjID,corrected_Brain_PAD = c(cPAD_rvm),Status=rep(c("Control","Patient"),c(nrow(Test_HC_data),nrow(Test_Patient_dat))))
dat$status_sch <- dat$status_trs <- dat$status_ssd <- rep(NA,nrow(dat))
dat$status_sch[test_dat$SubjID %in% Test_HC_data$SubjID ] <- "Control"
dat$status_sch[test_dat$SubjID %in% Test_schizophrenia_dat$SubjID ] <- "Schizophrenia"
dat$status_trs[test_dat$SubjID %in% Test_TRS_dat$SubjID ] <- "TRS"
dat$status_trs[test_dat$SubjID %in% Test_HC_data$SubjID ] <- "Control"
dat$status_ssd[test_dat$SubjID %in% Test_FESSD_dat$SubjID ] <- "FE-SSD"
dat$status_ssd[test_dat$SubjID %in% Test_HC_data$SubjID ] <- "Control"
dat %>% dplyr::select(corrected_Brain_PAD,Status) %>% group_by(Status) %>% summarise_each(funs(mean,sd,std.error))## # A tibble: 2 × 4
## Status mean sd std.error
## <chr> <dbl> <dbl> <dbl>
## 1 Control -2.09 6.47 0.447
## 2 Patient 1.87 8.76 0.574dat %>% dplyr::select(corrected_Brain_PAD,status_sch) %>% group_by(status_sch) %>% summarise_each(funs(mean,sd,std.error))## # A tibble: 3 × 4
## status_sch mean sd std.error
## <chr> <dbl> <dbl> <dbl>
## 1 Control -2.09 6.47 0.447
## 2 Schizophrenia 1.80 8.89 0.635
## 3 <NA> 2.25 8.15 1.34dat %>% dplyr::select(corrected_Brain_PAD,status_trs) %>% group_by(status_trs) %>% summarise_each(funs(mean,sd,std.error))## # A tibble: 3 × 4
## status_trs mean sd std.error
## <chr> <dbl> <dbl> <dbl>
## 1 Control -2.09 6.47 0.447
## 2 TRS 3.17 9.04 1.33
## 3 <NA> 1.55 8.68 0.635dat %>% dplyr::select(corrected_Brain_PAD,status_ssd) %>% group_by(status_ssd) %>% summarise_each(funs(mean,sd,std.error))## # A tibble: 3 × 4
## status_ssd mean sd std.error
## <chr> <dbl> <dbl> <dbl>
## 1 Control -2.09 6.47 0.447
## 2 FE-SSD 1.66 7.88 0.929
## 3 <NA> 1.97 9.14 0.721ggplot(dat,aes(x = corrected_Brain_PAD, fill = Status)) + geom_density(alpha = 0.5) + scale_fill_manual(values=palette) +
coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))
ggplot(dat[!is.na(dat$status_sch),],aes(x = corrected_Brain_PAD, fill = status_sch)) + scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))
ggplot(dat[!is.na(dat$status_ssd),],aes(x = corrected_Brain_PAD, fill = status_ssd)) + scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))
ggplot(dat[!is.na(dat$status_trs),],aes(x = corrected_Brain_PAD, fill = status_trs)) +scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))
library(ggplot2)
ggplot(dat,aes(x = corrected_Brain_PAD, fill = Status)) + geom_density(alpha = 0.5) + scale_fill_manual(values=palette[c(3,5)]) +
coord_cartesian(xlim =c(-30,30)) + theme_bw()
A <- ggplot(dat[!is.na(dat$status_sch),],aes(x = corrected_Brain_PAD, fill = status_sch)) + scale_fill_manual(values=palette[c(3,4)]) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))
B <- ggplot(dat[!is.na(dat$status_ssd),],aes(x = corrected_Brain_PAD, fill = status_ssd)) + scale_fill_manual(values=palette[c(3,1)]) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))
C <- ggplot(dat[!is.na(dat$status_trs),],aes(x = corrected_Brain_PAD, fill = status_trs)) +scale_fill_manual(values=palette[c(3,2)]) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))
library(ggpubr)
png("S6_2(RVR).png", width = 1400, height = 300)
ggarrange(A,B,C,
labels = c("A","B","C"),
ncol = 3, nrow = 1,
font.label = list(size=20))
dev.off()## png
## 2# ANOVA results (Control vs FE-SSD vs TRS)
dat_for_regression <- dat
replicated <- setdiff(which(!is.na(dat$status_trs) & !is.na(dat$status_ssd)),1:214)
dat[replicated,]## SubjID corrected_Brain_PAD Status status_ssd status_trs status_sch
## 415 T19 3.254155 Patient FE-SSD TRS Schizophrenia
## 419 T24 -1.427691 Patient FE-SSD TRS Schizophreniadat_for_regression <- rbind(dat_for_regression,dat[replicated,])
dat_for_regression[replicated,]$status_ssd <- NA
dat_for_regression[(nrow(dat_for_regression)-1):nrow(dat_for_regression),]$status_trs <- NA
table(dat_for_regression$status_trs)##
## Control TRS
## 209 46table(dat_for_regression$status_ssd)##
## Control FE-SSD
## 209 72test_dat_for_regression <- test_dat
test_dat_for_regression <- rbind(test_dat_for_regression,test_dat[replicated,])
dat_for_regression$status_ssd[is.na(dat_for_regression$status_ssd)] = dat_for_regression$status_trs[is.na(dat_for_regression$status_ssd)]
summary(res_aov <- aov(dat_for_regression$corrected_Brain_PAD~
dat_for_regression$status_ssd))## Df Sum Sq Mean Sq F value Pr(>F)
## dat_for_regression$status_ssd 2 1483 741.3 14.31 1.11e-06 ***
## Residuals 324 16786 51.8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 결측으로 인하여 117개의 관측치가 삭제되었습니다.pairwise.t.test(dat_for_regression$corrected_Brain_PAD,
dat_for_regression$status_ssd, p.adj = "fdr")$p.value## Control FE-SSD
## FE-SSD 2.479033e-04 NA
## TRS 3.036364e-05 0.2675513pairwise.t.test(dat_for_regression$corrected_Brain_PAD,
dat_for_regression$status_ssd, p.adj = "bonferroni")$p.value## Control FE-SSD
## FE-SSD 4.958066e-04 NA
## TRS 3.036364e-05 0.802654# t.test results (Control vs Schizo)
t.test(dat$corrected_Brain_PAD~dat$status_sch)##
## Welch Two Sample t-test
##
## data: dat$corrected_Brain_PAD by dat$status_sch
## t = -5.0093, df = 354.77, p-value = 8.63e-07
## alternative hypothesis: true difference in means between group Control and group Schizophrenia is not equal to 0
## 95 percent confidence interval:
## -5.416722 -2.362566
## sample estimates:
## mean in group Control mean in group Schizophrenia
## -2.088031 1.801613t.test(dat$corrected_Brain_PAD~dat$status_sch, var.equal=T)##
## Two Sample t-test
##
## data: dat$corrected_Brain_PAD by dat$status_sch
## t = -5.0592, df = 403, p-value = 6.411e-07
## alternative hypothesis: true difference in means between group Control and group Schizophrenia is not equal to 0
## 95 percent confidence interval:
## -5.401066 -2.378221
## sample estimates:
## mean in group Control mean in group Schizophrenia
## -2.088031 1.801613summary(res_aov <- aov(dat$corrected_Brain_PAD~dat$status_sch))## Df Sum Sq Mean Sq F value Pr(>F)
## dat$status_sch 1 1530 1530.3 25.59 6.41e-07 ***
## Residuals 403 24094 59.8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 결측으로 인하여 37개의 관측치가 삭제되었습니다.